Note. Boxplots display the interquartile range (IQR, center box), and the whiskers extend 1.5*IQR from the lower and upper hinge. The white point indicates the mean and the white center line indicates the median.


Data Preparation

In an initial preparatory step, we import the data into the R project environment and prepare the variables for further processing and later analyses.

Data Import

The data were collected using two different survey tools. For the study with sojourners (Study 1: worker) we used the survey platform Qualtrics XM, whereas the studies with international students (Study 2: student), and the international medical professionals (Study 3: medical) were conducted using the survey framework FormR. This means that the datasets had inconsistent file formats and naming conventions. For the Qualtrics study we pre-processed some variables to ease the import process (for the syntax files see the SPS files in ‘data/S1_Workers/processed/cleaned’ and for the raw data files see ‘data/S1_Workers/raw’). For the two other studies, we import the raw csv files from their respective folders.

# workers
# initial data cleaning was done in SPSS (syntax files are available in "")
dtWorker <- list(
  raw.pre = read_spss("data/S1_Workers/processed/cleaned/MT - Pre-Measure - 06-15-2018.sav"),
  raw.post = read_spss("data/S1_Workers/processed/cleaned/MT - Post-Measure - 06-15-2018.sav"),
  raw.morning = read_spss("data/S1_Workers/processed/cleaned/MT - Morning - 06-15-2018.sav"),
  raw.afternoon = read_spss("data/S1_Workers/processed/cleaned/MT - Afternoon - 06-15-2018.sav")
)

# students
dtStudents <- list(
  raw.pre = read.csv(file = "data/S2_Students/raw/AOTS_Pre.csv", header = T, sep = ","),
  raw.post = read.csv(file = "data/S2_Students/raw/AOTS_Post.csv", header = T, sep = ","),
  raw.daily = read.csv(file = "data/S2_Students/raw/AOTS_Daily.csv", header = T, sep = ",")
)

# young medical professionals
dtMedical <- list(
  raw.pre = c(),
  raw.post = c(),
  raw.daily = c()
)

Data Cleaning & Data Exclusions

Worker

For the sojourner sample data was collected in four separate surveys: (1) the pre-measurement, (2) the daily morning survey, (3) the daily afternoon survey, as well as (4) a post-measurement. We combine the four individual surveys into one cohesive dataframe and drop superfluous variables that are not relevant to the analyses relevant here. We then format the time and date variables and add person- and measurement indices (for easy and meaningful addressing of the data). We also exclude our own test data.
Note: All data preparation steps are saved in the ‘dtWorker’ list.

#  important names for Morning and Afternoon
names.m <- c(
  "StartDate",
  "EndDate",
  "Finished",
  "Duration__in_seconds_",
  "RecordedDate",
  "ExternalReference",
  "Meta_Operating_System",
  "Contact_dum",
  "number",
  "time",
  "duration_1",
  "dyad.group",
  "gr_size",
  "gr_type_1",
  "gr_type_2",
  "gr_type_3",
  "gr_type_4",
  "gr_type_5",
  "gr_type_6",
  "gr_type_7",
  "gr_type_8",
  "gr_type_9",
  "gr_type_10",
  "gr_type_11",
  "gr_type_12",
  "gr_type_13",
  "gr_type_14",
  "gr_type_15",
  "gr_type_16",
  "gr_type_17_TEXT",
  "gr_context_1",
  "gr_context_2",
  "gr_context_3",
  "gr_context_4",
  "gr_context_5",
  "gr_context_6",
  "gr_context_7",
  "gr_context_8",
  "gr_context_9",
  "gr_context_10",
  "gr_context_11",
  "gr_context_12",
  "gr_context_13_TEXT",
  "gr_context_14_TEXT",
  "gr_dutchness",
  "dyad_type_1",
  "dyad_type_2",
  "dyad_type_3",
  "dyad_type_4",
  "dyad_type_5",
  "dyad_type_6",
  "dyad_type_7",
  "dyad_type_8",
  "dyad_type_9",
  "dyad_type_10",
  "dyad_type_11",
  "dyad_type_12",
  "dyad_type_13",
  "dyad_type_14",
  "dyad_type_15",
  "dyad_type_16",
  "dyad_type_17_TEXT",
  "Context_1",
  "Context_2",
  "Context_3",
  "Context_4",
  "Context_5",
  "Context_6",
  "Context_7",
  "Context_8",
  "Context_9",
  "Context_10",
  "Context_11",
  "Context_12",
  "Context_13_TEXT",
  "Context_14_TEXT",
  "keyMotive",
  "keymotive_fulfillemt_1",
  "keyMotive_Dutch_1",
  "autonomy_1",
  "competence_1",
  "relatedness_self_1",
  "relatedness_other_1",
  "qualityAccidental_1",
  "qualityVoluntary_1",
  "qualityCooperative_1",
  "qualityDutchy_1",
  "quality_overall_1",
  "quality_meaning_1",
  "quality_star_1",
  "wantInt",
  "desire_type_1",
  "desire_type_2",
  "desire_type_3",
  "desire_type_4",
  "desire_type_5",
  "desire_type_6",
  "desire_type_7",
  "desire_type_8",
  "desire_type_9",
  "desire_type_10",
  "desire_type_11",
  "desire_type_12",
  "desire_type_13",
  "desire_type_14",
  "desire_type_15",
  "desire_type_16",
  "desire_type_17_TEXT",
  "desire_context_1",
  "desire_context_2",
  "desire_context_3",
  "desire_context_4",
  "desire_context_5",
  "desire_context_6",
  "desire_context_7",
  "desire_context_8",
  "desire_context_9",
  "desire_context_10",
  "desire_context_11",
  "desire_context_12",
  "desire_context_13_TEXT",
  "desire_context_14_TEXT",
  "Reason_nodesire",
  "keyMotive_noInt",
  "keyMotive_noInt_fulf_1",
  "autonomy_NoInt_1",
  "competence_NoInt_1",
  "relatedness_1_NoInt_1",
  "thermometerDutch_1",
  "thermometerDutchInt_2",
  "ExWB_1",
  "alertness1",
  "calmness1",
  "valence1",
  "alertness2",
  "calmness2",
  "valence2",
  "inNonDutch",
  "NonDutchNum",
  "NonDutchType_1",
  "NonDutchType_2",
  "NonDutchType_3",
  "NonDutchType_4",
  "NonDutchType_5",
  "NonDutchType_6",
  "NonDutchType_7",
  "NonDutchType_8",
  "NonDutchType_9",
  "NonDutchType_10",
  "NonDutchType_11",
  "NonDutchType_12",
  "NonDutchType_13",
  "NonDutchType_14",
  "NonDutchType_15_TEXT",
  "date",
  "time.0",
  "LocationLatitude",
  "LocationLongitude"
)

names.a <- c(names.m, "keyInteraction_1", "keyInteractionTime")

# Create reduced data sets for morning and afternoon
dat.mo <- dtWorker$raw.morning[, names.m]
dat.mo$daytime <- "morning"

dat.af <- dtWorker$raw.afternoon[, names.a]
dat.af$daytime <- "afternoon"

# merge morning and afternoon measurements with indicator [+ clean up]
daily.dat <- rbind.fill(dat.mo, dat.af)
daily.dat <- daily.dat[daily.dat$ExternalReference != 55951, ]
dtWorker$daily <- daily.dat
rm(dat.mo, dat.af, names.m, names.a, daily.dat)


# names for pre-measurement
names.pre <- c(
  "Finished",
  "age",
  "Gender",
  "Living",
  "roommate_1",
  "roommate_2",
  "roommate_3",
  "nationality",
  "SecondNationality",
  "timeNL_1",
  "Reason_2",
  "Reason_5",
  "Reason_7",
  "Reason_8_TEXT",
  "DutchLang",
  "occupation_1",
  "occupation_2",
  "occupation_3",
  "occupation_4",
  "occupation_7",
  "CurrentEducation_1",
  "education_level",
  "EduLang_2",
  "RUG_faculty",
  "Study.0",
  "association",
  "DutchMeetNum",
  "DutchFriends_1",
  "assimilation",
  "separation",
  "integration",
  "marginalization",
  "VIA_heritage",
  "VIA_Dutch",
  "SSAS_surrounding",
  "SSAS_privat",
  "SSAS_public",
  "autonomy",
  "relatedness",
  "competence",
  "anxiety",
  "swl",
  "alertness",
  "calmness",
  "valence",
  "date",
  "time",
  "City",
  "ZIP",
  "id"
)

# reduced data set for pre measurement
dat.pre.red <- dtWorker$raw.pre[, names.pre]

# merge with daily data [+ clean up]
df.pre <- merge(
  x = dtWorker$daily,
  y = dat.pre.red,
  by.x = "ExternalReference",
  by.y = "id",
  all = T
)
rm(names.pre)

# adjust duplicate names to fit to indicate daily or pre measurement
names(df.pre) <- gsub("[[:punct:]]x", ".daily", names(df.pre))
names(df.pre) <- gsub("[[:punct:]]y", ".pre", names(df.pre))

# names for post measurement
names.post <- c(
  "ExternalReference",
  "assimilation",
  "separation",
  "integration",
  "marginalization",
  "VIA_heritage",
  "VIA_Dutch",
  "anxiety",
  "swl",
  "rosenberg",
  "social_support",
  "stress",
  "discrimination",
  "discrimination_month",
  "NLE_1month",
  "NLE_6month",
  "NLE_12month"
)

# reduced data set for post-measurement
dat.post.red <- dtWorker$raw.post[, names.post]

# merge post measurement with pre- and daily data
df <- merge(
  x = df.pre,
  y = dat.post.red,
  by.x = "ExternalReference",
  by.y = "ExternalReference",
  all = T
)

# adjust duplicate names to indicate pre or post
names(df) <- gsub("[[:punct:]]x", ".pre", names(df))
names(df) <- gsub("[[:punct:]]y", ".post", names(df))

# add to list
dtWorker$combined <- df

# create data frame with cleaned data
df <- dtWorker$combined %>%
  filter(
    Finished.pre == 1,
    Finished.daily == 1,
    !is.na(ExternalReference)
  )

# add running number as measurement ID within participants
df$measureID <- rowidv(df, cols = c("ExternalReference"))

df <- df %>%
  mutate(
    PID = as.numeric(factor(ExternalReference)),
    # participant ID
    TID = measureID - 1,
    # time ID with t0 = 0 for meaningfull intercept interpretations
    date = substr(StartDate, 1, 10),
    # awkward way of extracting date (best converted to )
    time = substr(StartDate, 12, 19),
    # awkward way of extracting time
    daynum = as.numeric(factor(date)),
    # all days as numeric for ordering
    daycor = ifelse(
      daytime == "morning" &
        period_to_seconds(hms(time)) < period_to_seconds(hms("12:00:00")) |
        daytime == "afternoon" &
          period_to_seconds(hms(time)) < period_to_seconds(hms("19:00:00")),
      daynum - 1,
      daynum
    ),
    # correctly identify which date the questionnaire is about
    daycor.lead = sprintf("%02d", daycor),
    daytime.lt = ifelse(daytime == "morning", "a", "b"),
    # morning / afternoon to a / b
    day_time = paste(daycor.lead, daytime.lt, sep = "_"),
    # combine day id with morning / afternoon
    session = as.numeric(factor(day_time)),
    # day and time identifier as numeric id
    SubTime = chron::times(time.0),
    time.daily = as.character(time.daily),
    PPDate = as.Date(df$date.daily),
    number = replace_na(number, 0),
    NonDutchNum = replace_na(NonDutchNum, 0)
  )

dtWorker$clean <- df

# clean up
rm(df.pre, names.post, dat.post.red, dat.pre.red, df)

# Export reduced Data
# write.csv(dtWorker$clean, "data/processed/MT_clean-merged_07-05-2018.csv", row.names = F)
# save(dtWorker$clean, file = "data/processed/MT_clean-merged_07-05-2018.RData")

Student

For the student sample data was, similarly, collected in three separate surveys: (1) the pre-measurement, (2) the daily survey sent out at lunch and dinner time, and (3) a post-measurement. We combine the three individual surveys into one large dataframe and drop superfluous variables that are not relevant to the analyses relevant here. We exclude our own test data as well as one participant who entered the study twice (but gave different responses during the pre-measurement). We also reformat missing values and format core ID variables.
Note: All data preparation steps are saved in the ‘dtStudents’ list.

# our own test IDs
ownIDs <- c(
  "beautifulLionfishXXXR5rcgVBzGu8hPvOqrK8UBJBw4owvi9nfRFSFu3lMzYhE",
  "niceDogoXXXmB8JI5SFu78SF3DVof84mGUPPNUr14p2HYFTtp31a6D1OwAzM6F-K",
  "amusedQuailXXXmhuc_fpTp8vPkMwDH1BzjaH1d1kHSO1bsPEfsnaEYk4WeVBfPi",
  "juwGAbtXX0_1kmZtSVqKh3PGaHOICqUyU4iBkrT3nDsI_uifuD1gzKcZerxaM5FL"
)

# Prepare dfs for Cleaning
df.pre <- dtStudents$raw.pre %>%
  mutate_all(na_if, "") %>%
  mutate_all(na_if, "NA") %>%
  filter(!is.na(ended)) %>% # remove all who did not finish
  filter(!e_mail %in% .$e_mail[duplicated(.$e_mail)]) %>% # remove all who did the pre questionnaire multiple times (b/c inconsistent ratings scales)
  filter(!session %in% ownIDs) %>% # remove our own test
  mutate(session = as.character(session)) # turn factor into character strings (probably just precaution)

df.post <- dtStudents$raw.post %>%
  mutate_all(na_if, "") %>%
  mutate_all(na_if, "NA") %>%
  filter(!is.na(session)) %>% # remove own test runs
  filter(!session %in% ownIDs) %>% # remove our own test
  filter(session %in% df.pre$session) %>% # remove anyone who wasn't in the pre
  filter(!is.na(ended)) %>% # remove all who never finished
  filter(!session %in% .$session[duplicated(.$session)]) %>% # remove all duplicate sessions
  mutate(session = as.character(session)) # turn factor into character strings (probably just precaution)

df.daily <- dtStudents$raw.daily %>%
  mutate_all(na_if, "") %>%
  mutate_all(na_if, "NA") %>%
  filter(!session %in% ownIDs) %>% # remove our own test
  filter(session %in% df.pre$session) %>% # remove anyone who wasn't in the pre
  filter(!is.na(ended)) %>% # remove all who never finished
  mutate(session = as.character(session)) # turn factor into character strings (probably just precaution)

# merge daily with pre
dfPreDaily <- merge(
  x = df.daily,
  y = df.pre,
  by = "session",
  suffixes = c(".daily", ".pre"),
  all = F
)

# merge daily with post
dfCombined <- merge(
  x = dfPreDaily,
  y = df.post,
  by = "session",
  suffixes = c(".pre", ".post"),
  all = F
)

# add to list
dtStudents$clean <- dfCombined

# clean up workspace
rm(df.pre, df.daily, df.post, dfPreDaily, dfCombined, ownIDs)

Medical

To be done after preregistration.

# TBD

Calculate needed transformations

Worker

For the worker sample, the data transformation stage had three main aims:

  1. We first corrected time indicators within the surveys. In some cases participants completed their daily diary surveys for the afternoon after midnight. In these cases the measurement still is in reference to the previous day and is indicated in the corrected variable.
  2. We then created indices of scales. Some indices were multi-item scales while some indices combine equivalent measurement for different situational circumstances (e.g., competence perceptions after interactions and at measurement occasions without interactions).
  3. Finally, we calculated several basic participant summaries (averages across all measurement occasions).
df <- dtWorker$clean

# Time and Date Variables
# remove seconds from afternoon time
df$SubTime[df$daytime == "afternoon"] <- paste0(substring(as.character(df$time.0[df$daytime == "afternoon"]), 4, 8), ":00")
df$time.daily[df$daytime == "afternoon" &
  !is.na(df$time.daily != "<NA>")] <- paste0(substring(as.character(df$time.daily[df$daytime == "afternoon" &
  !is.na(df$time.daily != "<NA>")]), 4, 8), ":00")

# Correct morning / afternoon date where survey was collected the day after to indicate the correct date that was targeted
df$PPDate[df$SubTime < "11:50:00" &
  df$daytime == "morning"] <- df$PPDate[df$SubTime < "11:50:00" &
  df$daytime == "morning"] - 1
df$PPDate[df$SubTime < "18:50:00" &
  df$daytime == "afternoon"] <- df$PPDate[df$SubTime < "18:50:00" &
  df$daytime == "afternoon"] - 1

# Need scales
df$keyMotiveFulfilled <- rowSums(df[, c("keymotive_fulfillemt_1", "keyMotive_noInt_fulf_1")], na.rm = T)
df$autonomy.daily.all <- rowSums(df[, c("autonomy_1", "autonomy_NoInt_1")], na.rm = T)
df$competence.daily.all <- rowSums(df[, c("competence_1", "competence_NoInt_1")], na.rm = T)
# cor(df$relatedness_other_1, df$relatedness_self_1,use="complete.obs")
df$relatedness.daily.all <- rowMeans(df[, c(
  "relatedness_other_1",
  "relatedness_self_1",
  "relatedness_1_NoInt_1"
)], na.rm = T)

pairs.panels.new(
  df[c("relatedness_self_1", "relatedness_other_1")],
  labels = c(
    "I shared information about myself.",
    "X shared information about themselves."
  )
)

df$relatedness_1 <- rowMeans(df[, c("relatedness_other_1", "relatedness_self_1")], na.rm = T)

# summarize by participant (check that everything is within pp might not be the case for )
between <- df %>%
  group_by(ExternalReference) %>%
  mutate(
    CtContactNL = sum(Contact_dum),
    CtContactNonNl = sum(inNonDutch),
    CtContactNLAll = sum(number),
    CtContactNonNlAll = sum(NonDutchNum),
    AvKeyNeed = mean(keyMotiveFulfilled, na.rm = T),
    AvKeyNeedInt = mean(keymotive_fulfillemt_1, na.rm = T),
    AvKeyNeedNoInt = mean(keyMotive_noInt_fulf_1, na.rm = T),
    AvAutonomy = mean(autonomy.daily.all, na.rm = T),
    AvCompetence = mean(competence.daily.all, na.rm = T),
    AvRelatedness = mean(relatedness.daily.all, na.rm = T),
    AvThermo = mean(thermometerDutch_1, na.rm = T),
    AvWB = mean(ExWB_1, na.rm = T)
  ) %>%
  ungroup() %>%
  mutate(
    CtContactNL_c = scale(CtContactNL, scale = FALSE),
    AvKeyNeedInt_c = scale(AvKeyNeedInt, scale = FALSE),
    AvKeyNeed_c = scale(AvKeyNeed, scale = FALSE),
    CtContactNL_z = scale(CtContactNL, scale = TRUE),
    AvKeyNeedInt_z = scale(AvKeyNeedInt, scale = TRUE),
    AvKeyNeed_z = scale(AvKeyNeed, scale = TRUE)
  )

warning(
  "some variable transformations (esp. _c and _z) might be across all participants (i.e., not within PP)."
)

dtWorker$full <- between
rm(df, between)

# Between participants contact frequency
workerContactFreq <- dtWorker$full %>%
  group_by(ExternalReference) %>%
  summarise(
    n = n(),
    SumContactNL = sum(Contact_dum),
    PercContactNL = SumContactNL / n * 100,
    SumContactNLAll = sum(number),
    AvAttitude = mean(thermometerDutch_1, na.rm = T)
  ) %>%
  mutate(
    WinSumContactNL = DescTools::Winsorize(SumContactNL),
    WinSumContactNLAll = DescTools::Winsorize(SumContactNLAll)
  )

# dataframe where interaction types are recoded
workerInteractionType <- dtWorker$full %>%
  mutate(
    OutgroupInteraction = as_factor(Contact_dum),
    NonOutgroupInteraction = as_factor(inNonDutch)
  )

# save cleaned data
# save(df.btw, file = "data/processed/df.btw.RData")
# write_sav(df.btw, "data/processed/MT_clean-merged_pre-post.sav")

# export data to Mplus
# df.mplus = remove_all_labels(select(df,
#                                     PID, session,
#                                     thermometerDutch_1, inNonDutch, Contact_dum,
#                                     keyMotiveFulfilled, autonomy.daily.all, competence.daily.all, relatedness.daily.all))
# names(df.mplus)= c("PID", "session", "att", "intin", "intout", "keymot", "aut", "comp", "rel")
# mplus = df.mplus[order(df.mplus$PID, df.mplus$session),]
# mplus.intcont = mplus[mplus$intout==1,]
# prepareMplusData(mplus.intcont, "data/processed/dynamic-subset-intonly.dat")

Student

For the worker sample, the data transformation stage had three main aims:

  1. We first create person, survey type, and measurement ID variables.
  2. We then created indices of scales. Some indices were multi-item scales while some indices combine equivalent measurement for different situational circumstances (e.g., competence perceptions after interactions and at measurement occasions without interactions).
  3. We add information about the interaction partner to the beep during which a person was selected as an interaction partner.
  4. We cluster mean-center key variables within participants.
  5. Finally, we calculated several basic participant summaries (averages across all measurement occasions).
df <- dtStudents$clean

# Add ID variables
df$PID <- as.numeric(factor(df$session)) # participant ID

# order time
df$TID <-
  factor(df$date_period, levels = unique(dtStudents$raw.daily$date_period))
df$TIDnum <- as.numeric(df$TID) # get numeric TID

# check whether time ordering worked
df <- df %>%
  arrange(PID, TID) # %>%
# View()

# Interaction as Factor
df$interaction.f <-
  factor(df$Interaction,
    levels = c("no interaction", "Dutch", "Non-Dutch")
  )
df$intNL <- ifelse(df$Interaction == "Dutch", 1, 0)
df$intNonNL <- ifelse(df$Interaction == "Non-Dutch", 1, 0)

# -------------------------------------------------------------------------------------------------------------
#                                       Combine Variables
# -------------------------------------------------------------------------------------------------------------
# Relatedness
pairs.panels.new(
  df[c("RelatednessSelf", "RelatednessOther")],
  labels = c(
    "I shared information about myself.",
    "X shared information about themselves."
  )
)

df$RelatednessInteraction <-
  rowMeans(df[c("RelatednessSelf", "RelatednessOther")], na.rm = T)
df$RelatednessInteraction[df$RelatednessInteraction == "NaN"] <-
  NA
# Relatedness Overall (JANNIS NOT SURE THESE ARE CORRECT, CHANGE ROWS?; J: Changed "NaN" in df$RelatednessInteraction to NA() should work now)
df$Relatedness <-
  rowMeans(df[, c("RelatednessInteraction", "RelatednessNoInteraction")],
    na.rm =
      T
  )
# Pro-Sociality
df$ProSo <-
  rowMeans(df[, c("ProSo1", "ProSo2", "ProSo3", "ProSo4")], na.rm = T)
# Anti-Sociality
df$AntiSo <-
  rowMeans(df[, c("AntiSo1", "AntiSo2", "AntiSo3", "AntiSo4")], na.rm = T)


# -------------------------------------------------------------------------------------------------------------
#                                 Add Variables related to interaction partner
# -------------------------------------------------------------------------------------------------------------
# create function for later lapply
createIntPartDf <- function(inp) {
  # prepare the dataframe so that we can forloop over it later
  tmp <- data.frame(
    CC = as.character(inp$CC),
    NewCC = as.character(inp$NewCC),
    NewName = as.character(inp$NewName),
    NewCloseness = inp$NewCloseness,
    NewGender = inp$NewGender,
    NewEthnicity = as.character(inp$NewEthnicity),
    NewRelationship = as.character(inp$NewRelationship)
  )

  tmp$CC2 <- recode(tmp$CC, "SOMEONE ELSE" = "NA")
  tmp$CC2 <-
    ifelse(
      tmp$CC == 1 |
        tmp$CC == "SOMEONE ELSE",
      as.character(tmp$NewName),
      as.character(tmp$CC2)
    )
  # maybe add [[:space:]]\b to remove space before word boundary or ^[[:space:]] to remove space in the beginning of a string
  tmp$CC2 <- gsub("^[[:space:]]", "", tmp$CC2)
  tmp$NewName <- gsub("^[[:space:]]", "", tmp$NewName)

  # open the variables that will be filled up in the foor-loop
  tmp$closeness <- rep(NA, nrow(tmp))
  tmp$gender <- rep(NA, nrow(tmp))
  tmp$ethnicity <- rep(NA, nrow(tmp))
  tmp$relationship <- rep(NA, nrow(tmp))

  # Run the for-loop. It finds the variables related to the name of the interaction partner. If there is a repeating interaction
  # partner (i.e. CC2) it takes the value (i.e. NewCloseness) from the first interaction (i.e. NewName)
  for (i in 1:nrow(tmp)) {
    if (is.na(tmp$CC2[i])) {
      next
    } else {
      tmp$closeness[i] <-
        na.omit(tmp$NewCloseness[as.character(tmp$CC2[i]) == as.character(tmp$NewName)])[1] # find closeness where CC2 matches NewName (na.omit + [1] to get the number)
      tmp$gender[i] <-
        na.omit(tmp$NewGender[as.character(tmp$CC2[i]) == as.character(tmp$NewName)])[1] # (na.omit + [1] to get the number and not the rest of the na.omit list)
      tmp$ethnicity[i] <-
        na.omit(as.character(tmp$NewEthnicity[as.character(tmp$CC2[i]) == as.character(tmp$NewName)]))[1] # PROBLEM IS THAT THERE ARE TOO MANY NA's: Difficult to deal with
      tmp$relationship[i] <-
        na.omit(as.character(tmp$NewRelationship[as.character(tmp$CC2[i]) == as.character(tmp$NewName)]))[1]
    }
  }

  out <- tmp
  out
}

# split df per participants and run function
PP <- split(df, df$PID)
PP <- lapply(PP, createIntPartDf)
rm(createIntPartDf)

# add variables back to df
remergePP <- do.call(rbind.data.frame, PP)
colnames(remergePP) <-
  paste(colnames(remergePP), "_Calc", sep = "")
df <- cbind(df, remergePP)
rm(remergePP, PP)

# -------------------------------------------------------------------------------------------------------------
#                                 Center Relevant Variables
# -------------------------------------------------------------------------------------------------------------

df <- df %>%
  group_by(PID) %>%
  mutate(
    KeyNeedFullfillment.cm = mean(KeyNeedFullfillment, na.rm = TRUE),
    # cluster mean (mean of PP)
    KeyNeedFullfillment.cwc = KeyNeedFullfillment - KeyNeedFullfillment.cm,
    # cluster mean centered (within PP centered)
    closeness.cm = mean(closeness_Calc, na.rm = TRUE),
    closeness.cwc = closeness_Calc - closeness.cm
  ) %>%
  ungroup()

# store
dtStudents$full <- df
rm(df)

# Between participants contact frequency
studentContactFreq <- dtStudents$full %>%
  group_by(PID) %>%
  summarise(
    n = n(),
    SumContactNL = sum(InteractionDumDutch),
    PercContactNL = SumContactNL / n * 100,
    SumContactNLAll = sum(ContactNum[InteractionDumDutch == 1], na.rm = TRUE),
    AvAttitude = mean(AttitudesDutch, na.rm = TRUE),
    AvQuality = mean(quality_overall, na.rm = TRUE)
  ) %>%
  mutate(
    WinSumContactNL = DescTools::Winsorize(SumContactNL),
    WinSumContactNLAll = DescTools::Winsorize(SumContactNLAll)
  )

# dataframe where interaction types are recoded
studentInteractionType <- dtStudents$full %>%
  mutate(
    NonDutchContact = replace_na(NonDutchContact, 2), # make second non-Dutch countable
    NonDutchContact = NonDutchContact*-1+2 # recode (yes = 1 -> 1, no = 2 -> 0)
  ) %>%
  mutate(
    OutgroupInteraction = factor(
      InteractionDumDutch,
      levels = c(0, 1),
      labels = c("No", "Yes")
    ),
    NonOutgroupInteraction = factor(
      rowSums(select(., c(InteractionDumNonDutch, NonDutchContact))), # combine the two non-Dutch Q.,
      levels = c(0, 1),
      labels = c("No", "Yes")
    )
  )

# select a subset of IDs to display in plots
studentPltIDs <-
  studentInteractionType %>%
  group_by(PID) %>%
  summarise(n = n()) %>%
  slice_max(n, n = 20) %>% # chose the 20 with the most number of measurements
  select(PID) %>%
  as.matrix %>%
  as.vector

# select a subset of IDs to display in plots (only outgroup interactions)
studentOutPltIDs <-
  studentInteractionType %>%
  filter(OutgroupInteraction == "Yes") %>%
  group_by(PID) %>%
  summarise(n = n()) %>%
  slice_max(n, n = 20) %>% # chose the 20 with the most number of measurements
  select(PID) %>%
  as.matrix %>%
  as.vector

Medical

To be done after preregistration.

# TBD

Worker Sample

The first sample we assess is the smaller sojourner study. We will sequentially test our main hypotheses for this study:

  1. Based on the most general understanding of the contact hypothesis, an increase in frequency and quality of contact should jointly account for changes in more favorable outgroup attitudes.
    1. Participants with more intergroup interactions should have a more favorable outgroup attitudes.
    2. Outgroup attitudes should be higher after an intergroup interaction compared to a non-outgroup interaction.
    3. Participants with more intergroup interactions should have a more favorable outgroup attitudes depending on the average interaction quality.
  2. Based on our proposal, intergroup interactions with higher situational core need fulfillment should predict more favorable outgroup attitudes due to more positive interaction quality perceptions.
    1. Outgroup attitudes should be more favorable after intergroup interactions with high key need fulfillment.
    2. Interaction Quality should be perceived as more positive after intergroup interactions with higher key need fulfillment.
    3. The variance explained in outgroup attitudes by key need fulfillment should to a large extend be assumed by interaction quality.
    4. The effect of key need fulfillment on outgroup attitudes should be specific to intergroup interactions and not be due to need fulfillment in general. Thus, the effect of key need fulfillment on outgroup attitudes should stronger for intergroup interact than for ingroup interactions.
    5. The effect of key need fulfillment on outgroup attitudes should be persist even when taking other fundamental psychological needs into account. Thus, the effect of key need fulfillment on outgroup attitudes should remain strong even after controlling for autonomy, competence, and relatedness fulfillment during the interaction (cf., self-determination theory).

Data Description

Still in ‘scr/workerDescriptive.R’. Needs to be merged with this document.

Contact Hypothesis

We test the most general contact hypothesis in two steps. First, we assess whether more intergroup interactions are related to to more positive outgroup attitudes. Second, we test whether a potential positive effect on outgroup attitudes depends on the interaction quality (jointly with the number of interactions).

Interaction Frequency and Outgroup Attitudes

To test the impact of the overall number of interactions, we hope to find that there is a significant relationship between the number of interactions a participant had and the average outgroup attitude.

\[\begin{equation} \tag{1} r_{ContactFrequency, OutgroupAttitudes} \neq 0 \end{equation}\]

# correlation panel
pairs.panels.new(
  workerContactFreq %>% select(SumContactNL, SumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  workerContactFreq %>% select(WinSumContactNL, WinSumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

We find that neither the number of interactions nor the number of measurement beeps with an interaction are significantly related with the average outgroup attitudes. This is to say that within our data, participants with more outgroup interactions did not have significantly more positive outgroup attitudes. This might be due to the aggregation within the participants or the small sample size of between participant data. Nonetheless, the aggregate data does not support the notion that simply having more interactions with an outgroup results in more positive outgroup attitudes.

Outgroup Attitudes by Interaction Type

In a next step we take into account that having an interaction with an outgroup member, does not happen in a social vacuum. Participants who indicated that they had an interaction with an outgroup member include measurement occasions during which someone either only had an interaction with an outgroup member as well as times during which a person also had interaction(s) with a non-Dutch person. Inversely, participants who indicated that they did not have an interaction with a Dutch person might either have had no interaction at all or had an interaction with a non-Dutch person. We, thus consider all possible combinations and their independent influences on outgroup attitudes.

OLS regression

We first assess the impact of the different interaction types across all measurement points (lumping all beeps together).

\[\begin{equation} \tag{2} Attitude = OutgroupInteraction + NonOutgroupInteraction \end{equation}\]

# between participants interaction type
workerContactType <- dtWorker$full %>%
  group_by(
    Contact_dum,
    inNonDutch
  ) %>%
  summarise(
    n = n(),
    AttitudeM = mean(thermometerDutch_1, na.rm = T),
    AttitudeSD = sd(thermometerDutch_1, na.rm = T),
    AttitudeSE = AttitudeSD / sqrt(n),
    AttitudeLwr = AttitudeM - 1.96 * AttitudeSE,
    AttitudeUpr = AttitudeM + 1.96 * AttitudeSE
  ) %>%
  ungroup() %>%
  mutate(InteractionType = paste(
    ifelse(Contact_dum == 1, "Out", "NoOut"),
    ifelse(inNonDutch == 1, "In", "NoIn"),
    sep = ", "
  ))

# plot bar chart
ggplot(
  workerContactType,
  aes(
    y = AttitudeM,
    x = as_factor(Contact_dum),
    fill = as_factor(inNonDutch)
  )
) +
  geom_bar(
    stat = "identity",
    color = "black",
    position = position_dodge()
  ) +
  geom_errorbar(aes(ymin = AttitudeM, ymax = AttitudeUpr),
    width = .2,
    position = position_dodge(.9)
  ) +
  labs(
    fill = "Non-Outgroup Interaction",
    x = "Outgroup Interaction",
    y = "Outgroup Attitude",
    title = "Outgroup Attitudes by Interaction Type [95% CI]"
  ) +
  scale_fill_grey(
    start = 0.2,
    end = 0.8
  ) +
  scale_y_continuous(
    limits = c(0, 100),
    breaks = c(0, 15, 30, 40, 50, 60, 70, 85, 100),
    labels = c(
      "Very cold or unfavorable feelings 0°",
      "Quite cold and unfavorable feelings 15°",
      "Fairly cold and unfavorable feelings 30°",
      "A bit cold and unfavorable feelings 40°",
      "No feeling at all 50°",
      "A bit warm and favorable feelings 60°",
      "Fairly warm and favorable feelings 70° ",
      "Quite warm and favorable feelings 85° ",
      "Very warm and favorable feelings 100° "
    )
  ) +
  theme_Publication()

# create list to store Worker models
mdlWorker <- list()

# regression
mdlWorker$lmAttInt <-
  lm(thermometerDutch_1 ~ OutgroupInteraction * NonOutgroupInteraction,
    data = workerInteractionType
  )
# summary(lmWorkerAttInteraction)

summ(
  mdlWorker$lmAttInt,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type OLS linear regression
F(3,1221) 4.867
R² 0.012
Adj. R² 0.009
Est. 2.5% 97.5% t val. p
(Intercept) 69.507 67.899 71.116 84.796 0.000
OutgroupInteraction 3.620 0.378 6.862 2.191 0.029
NonOutgroupInteraction -0.513 -2.593 1.566 -0.484 0.628
OutgroupInteraction:NonOutgroupInteraction -0.098 -4.021 3.826 -0.049 0.961
Standard errors: OLS; Continuous predictors are mean-centered.

We find that while controlling for interactions with non-Dutch people, outgroup attitudes were significantly higher when participants had an interaction with a Dutch person. The effect is relatively small (3.62 points on a 0–100 scale). More importantly, however, this analysis lumps all ESM beeps from every participants together and ignores that the data is nested within participants.

ML regression

We, thus, proceed with a multilevel analysis, which acknowledges that the measurements are nested within participants.

Unconditional model

We start by creating an empty random intercept model (i.e., let the outgroup attitude intercept be different between participants; unconditional model).

\[\begin{equation} \tag{3} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorker$lmerAttNull <-
  lme4::lmer(thermometerDutch_1 ~ 1 + (1 | PID),
    data = dtWorker$full
  ) # use optim if it does not converge

mdlWorker$lmeAttNull <-
  lme(
    thermometerDutch_1 ~ 1,
    random = ~ 1 | PID,
    data = dtWorker$full,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(lmerWorkerAttNull) #or with the lme function
summ(mdlWorker$lmerAttNull, digits = 3, center = TRUE)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8805.880
BIC 8821.213
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.698
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 71.338 2.695 26.466 22.053 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.797
Residual 8.425
Grouping Variables
Group # groups ICC
PID 23 0.698 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1~1 + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlWorker$varAttNull <- 
  VarCorr(mdlWorker$lmeAttNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlWorker$tauAttNull <- 
  as.numeric(mdlWorker$varAttNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlWorker$sigmaAttNull <- 
  as.numeric(mdlWorker$varAttNull[2])
# The ICC estimate (between/between+within) is:
mdlWorker$IccAttNull <-
  (as.numeric(mdlWorker$varAttNull[1]) / (as.numeric(mdlWorker$varAttNull[1]) + as.numeric(mdlWorker$varAttNull[2])))
mdlWorker$IccPercAttNull <-
  ((as.numeric(mdlWorker$varAttNull[1]) / (as.numeric(mdlWorker$varAttNull[1]) + as.numeric(mdlWorker$varAttNull[2])))) * 100

We then compare the random intercept model to a model without a random intercept (i.e., without levels at all).

# Create and save Model
mdlWorker$glsAttNull  <-
  gls(thermometerDutch_1 ~ 1,
      data = dtWorker$full,
      control = list(opt = "nlmimb"))

# calculate Deviances manually:
mdlWorker$DevianceGlsNull <- logLik(mdlWorker$glsAttNull) * -2
mdlWorker$DevianceMlNull <- logLik(mdlWorker$lmeAttNull) * -2

# Compare the two null models:
anova(mdlWorker$glsAttNull,
      mdlWorker$lmeAttNull) %>% 
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 1: Worker: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
glsAttNull 1 2 10133 10144 -5065
lmeAttNull 2 3 8806 8821 -4400 1 vs 2 1329.406 < .001

Comparing the two empty model, we find that there is indeed a significant amount of variance explained by including a random intercept.

To assess the variances within and between participants we look at the ICC \(\tau_{00}^2 / (\tau_{00}^2 + \sigma^2)\): The ratio of the between-cluster variance to the total variance is called the Intraclass Correlation. It tells you the proportion of the total variance in Y that is accounted for by the clustering. (In this case the clustering means clustering observations per participant).

We find that an estimated 69.76% of the variation in Feeling Thermometer scores is explained by between participant differences (clustering by PID). This is to say that 69.76% of the variance in any individual report of Attitudes towards the Dutch can be explained by the properties of the individual who provided the rating. And we find that including ‘participant’ as a predictor adds significantly to the model.

random intercept with predictors

To this random intercept model we now add the two types of interactions possible at each measurement point as contemporaneous predictors of outgroup attitudes. That is: We check whether within participants having an outgroup interaction (or a non-outgroup interaction) is associated with more positive outgroup attitudes at the same measurement point.

\[\begin{equation} \tag{4} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorker$lmeInterceptAttType <-
  lme(
    thermometerDutch_1 ~ OutgroupInteraction + NonOutgroupInteraction,
    random =  ~ 1 | PID,
    data = workerInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    thermometerDutch_1 ~ OutgroupInteraction + NonOutgroupInteraction + (1 | PID),
    data = workerInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8788.218
BIC 8813.771
Pseudo-R² (fixed effects) 0.006
Pseudo-R² (total) 0.703
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.343 65.002 75.683 25.814 22.897 0.000
OutgroupInteraction 2.477 1.364 3.589 4.365 1204.135 0.000
NonOutgroupInteraction 0.427 -0.683 1.538 0.754 1204.911 0.451
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.811
Residual 8.362
Grouping Variables
Group # groups ICC
PID 23 0.701 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeAttNull, 
      mdlWorker$lmeInterceptAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 2: Worker: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeAttNull 1 3 8806 8821 -4400
lmeInterceptAttType 2 5 8788 8814 -4389 1 vs 2 21.663 < .001 
# Save variances
mdlWorker$varInterceptAttType <- 
  lme4::VarCorr(mdlWorker$lmeInterceptAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwInterceptAttType <-
  1 - (as.numeric(mdlWorker$varInterceptAttType[1]) / as.numeric(mdlWorker$varAttNull[1]))
mdlWorker$varBtwPercInterceptAttType <-
  (1 - (as.numeric(mdlWorker$varInterceptAttType[1]) / as.numeric(mdlWorker$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinInterceptAttType <-
  1 - (as.numeric(mdlWorker$varInterceptAttType[2]) / as.numeric(mdlWorker$varAttNull[2]))
mdlWorker$varWithinPercInterceptAttType <-
  (1 - (as.numeric(mdlWorker$varInterceptAttType[2]) / as.numeric(mdlWorker$varAttNull[2]))) * 100

We find that a random intercept model with the two interaction types as predictors explains significantly more variance then an empty random intercept model. Looking at the individual coefficients, we find that having an outgroup interaction is indeed associated with significantly more positive outgroup attitudes, while having an interaction with a non-Dutch person does not significantly relate to more positive or negative attitudes towards the Dutch.

TL;DR: Interaction with Dutch is great predictor, interactions with non-Dutch is not.

random slope

In a next step, we check whether further letting the effect of the different interaction types vary between participants explains additional variance in outgroup attitudes (i.e., random slopes).

\[\begin{equation} \tag{5} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorker$lmeSlopesAttType <- lme(
  thermometerDutch_1 ~
    OutgroupInteraction + NonOutgroupInteraction,
  random = ~ 1 + OutgroupInteraction + NonOutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = workerInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorker$lmerSlopesAttType <- lmer(
    thermometerDutch_1 ~
      OutgroupInteraction + NonOutgroupInteraction +
      (1 + OutgroupInteraction + NonOutgroupInteraction | PID),
    data = workerInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8793.695
BIC 8844.802
Pseudo-R² (fixed effects) 0.006
Pseudo-R² (total) 0.710
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.371 65.001 75.741 25.682 21.938 0.000
OutgroupInteraction 2.572 1.079 4.065 3.376 19.291 0.003
NonOutgroupInteraction 0.434 -0.693 1.562 0.755 205.804 0.451
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.883
PID OutgroupInteraction 2.248
PID NonOutgroupInteraction 0.424
Residual 8.307
Grouping Variables
Group # groups ICC
PID 23 0.706 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeInterceptAttType, 
      mdlWorker$lmeSlopesAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 3: Worker: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeInterceptAttType 1 5 8788 8814 -4389
lmeSlopesAttType 2 10 8794 8845 -4387 1 vs 2 4.165 0.526 
# Save variances
mdlWorker$varSlopesAttType <- 
  lme4::VarCorr(mdlWorker$lmeSlopesAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwSlopesAttType <- 
  1 - (as.numeric(mdlWorker$varSlopesAttType[1]) / as.numeric(mdlWorker$varInterceptAttType[1]))
mdlWorker$varBtwPercSlopesAttType <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttType[1]) / as.numeric(mdlWorker$varInterceptAttType[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinSlopesAttType <- 
  1 - (as.numeric(mdlWorker$varSlopesAttType[2]) / as.numeric(mdlWorker$varInterceptAttType[2]))
mdlWorker$varWithinPercSlopesAttType <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttType[2]) / as.numeric(mdlWorker$varInterceptAttType[2]))) * 100

# Assumption Checks:
mdlWorker$diagSlopesAttType <-
  sjPlot::plot_model(mdlWorker$lmerSlopesAttType, type = "diag")
grid.arrange(
  mdlWorker$diagSlopesAttType[[1]],
  mdlWorker$diagSlopesAttType[[2]]$`PID`,
  mdlWorker$diagSlopesAttType[[3]],
  mdlWorker$diagSlopesAttType[[4]]
)

# Plot prediction model
mdlWorker$predSlopesAttType <- 
  workerInteractionType %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorker$lmeSlopesAttType,
                           workerInteractionType,
                           re.form = NA
                           )
         )

(
  mdlWorker$predPltSlopesAttType <-
    ggplot(data = mdlWorker$predSlopesAttType, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Worker_PredictionPlot_SlopesAttType.png",
  mdlWorker$predPltSlopesAttType,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does not add significantly beyond the random intercept model. This is unusual because this might indicate the the effect is very consistent across participants. However, this might also be the case due to a small number of participants, or other measurement issues.

TL;DR: Random slopes don’t add much for this super simple model.

Interaction Frequency and Interaction Quality

In a final step we check whether the effect outgroup interactions, in part, depends on the quality during the interactions. Because we can only assess interaction quality when there is an interaction, it is difficult to assess this with the interaction dummy as a within person predictor. Instead, we will use an aggregate measure of interaction quality and average interaction quality to consider the two predictors jointly.

\[\begin{equation} \tag{6} Attitude = ContactFreq \times AverageContactQual \end{equation}\]

# prepare data
workerAvFreQual <- dtWorker$full %>%
  group_by(ExternalReference) %>%
  summarise(
    n = n(),
    SumContactNL = sum(Contact_dum),
    PercContactNL = SumContactNL / n * 100,
    SumContactNLAll = sum(number),
    AvAttitude = mean(thermometerDutch_1, na.rm = TRUE),
    AvQuality = mean(quality_overall_1, na.rm = TRUE)
  ) %>%
  mutate(
    WinSumContactNL = DescTools::Winsorize(SumContactNL),
    WinSumContactNLAll = DescTools::Winsorize(SumContactNLAll)
  )

# correlation panel
pairs.panels.new(
  workerAvFreQual %>% select(SumContactNL, SumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  workerAvFreQual %>% select(WinSumContactNL, WinSumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

Within the data, we find a medium sized correlation between the participants’ Average Interaction Quality and their Average Outgroup Attitudes. Thus within our data participants with a higher quality outgroup interactions also held more positive attitudes towards that group. However, the frequency of intergroup interactions had no meaningful correlation with either the average interaction quality or average outgroup attitudes.

# regression
mdlWorker$lmAttFreqQualX <-
  lm(AvAttitude ~ SumContactNL * AvQuality, data = workerAvFreQual)
# summary(lmWorkerAttFreqQualX)

summ(
  mdlWorker$lmAttFreqQualX,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 21 (2 missing obs. deleted)
Dependent variable AvAttitude
Type OLS linear regression
F(3,17) 6.663
R² 0.540
Adj. R² 0.459
Est. 2.5% 97.5% t val. p
(Intercept) 70.930 66.519 75.341 33.927 0.000
SumContactNL 0.269 -0.150 0.688 1.354 0.193
AvQuality 0.765 0.288 1.242 3.385 0.004
SumContactNL:AvQuality -0.049 -0.085 -0.014 -2.954 0.009
Standard errors: OLS; Continuous predictors are mean-centered. 
interactions::interact_plot(
  mdlWorker$lmAttFreqQualX,
  pred = AvQuality,
  modx = SumContactNL,
  interval = TRUE,
  partial.residuals = TRUE
)

interactions::johnson_neyman(mdlWorker$lmAttFreqQualX,
                             pred = AvQuality,
                             modx = SumContactNL,
                             alpha = .05)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SumContactNL is OUTSIDE the interval [23.66, 75.42], the slope of AvQuality is p < .05.
## 
## Note: The range of observed values of SumContactNL is [2.00, 51.00]

We find that interaction quality is significantly related to higher outgroup attitudes (albeit with a small effect size). We also find that in our sample with an increasing number of interactions the positive effect of interaction quality becomes weaker. However, it should be noted that this is based on data aggregating all within participant nuances and is only the date of 21 people.

Need Fulfillment

The main proposal of our article is that the success of an outgroup interaction might be explained by whether or not the interaction fulfilled the person’s core situational need. This should, in turn, be due to a higher perceived interaction quality. We will this sequentially test whether the fulfillment of the core need during an interaction is (1) related to more positive outgroup attitudes, (2) higher perceived interaction quality, and (3) whether the variance explained by the core need is assumed by the perceived interaction quality if considered jointly.

Need fulfillment and Attitudes

In a first step we, thus, test the relationship between outgroup attitudes and the fulfillment of the core situational need during the interaction.

Unconditional model

We again start by creating an empty random intercept model (i.e., let the outgroup attitude intercept be different between participants; unconditional model). Note that this unconditional model differs from the empty model created earlier because for this set of analyses we will only consider the subsample of measurement points during which an outgroup interaction was reported. This is necessary because measurements of needs during the interaction and perceived interaction quality are only meaningful within an interaction context.

# see how large our outgroup interaction subset actually is
tbl_cross(
  workerInteractionType,
  row = OutgroupInteraction,
  col = NonOutgroupInteraction,
  percent = "row"
)
Characteristic Did you meet non-Dutch people this morning? (in person for at least 10 minutes) Total
no yes
Did you meet a Dutch person this morning? (In person interaction for at least 10 minutes)
No 337 (40%) 501 (60%) 838 (100%)
Yes 110 (28%) 277 (72%) 387 (100%)
Total 447 (36%) 778 (64%) 1,225 (100%)
 
# create outgroup interaction subset
workerOutgroupInteraction <- workerInteractionType %>%
  filter(OutgroupInteraction == "Yes")

# create empty list to organize models
mdlWorkerOut <- 
  list(
    Att = list(),
    Qlt = list()
  )

Note to self: For completeness, we should probably check this against an empty model without random intercept.

\[\begin{equation} \tag{7} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmerNull <-
  lme4::lmer(thermometerDutch_1 ~ 1 + (1 | PID), 
             data = workerOutgroupInteraction) # use optim if it does not converge
mdlWorkerOut$Att$lmeNull <-
  lme(
    thermometerDutch_1 ~ 1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.ML.r) #or with the lme function
summ(mdlWorkerOut$Att$lmerNull, digits = 3, center = TRUE)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2863.460
BIC 2875.336
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.684
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 72.550 2.910 24.933 19.198 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.069
Residual 8.883
Grouping Variables
Group # groups ICC
PID 21 0.684 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1~1 + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlWorkerOut$Att$varNull <- 
  VarCorr(mdlWorkerOut$Att$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlWorkerOut$Att$tauNull <- 
  as.numeric(mdlWorkerOut$Att$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlWorkerOut$Att$sigmaNull <- 
  as.numeric(mdlWorkerOut$Att$varNull[2])
# The ICC estimate (between/between+within) is:
mdlWorkerOut$Att$IccNull <-
  (as.numeric(mdlWorkerOut$Att$varNull[1]) / (as.numeric(mdlWorkerOut$Att$varNull[1]) + as.numeric(mdlWorkerOut$Att$varNull[2])))
mdlWorkerOut$Att$IccPercNull <-
  ((as.numeric(mdlWorkerOut$Att$varNull[1]) / (as.numeric(mdlWorkerOut$Att$varNull[1]) + as.numeric(mdlWorkerOut$Att$varNull[2])))) * 100

random intercept with level one predictors

We thus add the core interaction need fulfillment to the multilevel random intercept model.

\[\begin{equation} \tag{8} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmeInterceptCore <-
  lme(
    thermometerDutch_1 ~ keymotive_fulfillemt_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(thermometerDutch_1 ~ keymotive_fulfillemt_1 + (1 | PID), 
       data = workerOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2842.675
BIC 2858.508
Pseudo-R² (fixed effects) 0.033
Pseudo-R² (total) 0.713
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 72.690 66.947 78.433 24.808 19.259 0.000
keymotive_fulfillemt_1 0.147 0.094 0.201 5.404 372.954 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.182
Residual 8.556
Grouping Variables
Group # groups ICC
PID 21 0.704 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Att$lmeNull, 
      mdlWorkerOut$Att$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 4: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeNull 1 3 2863 2875 -1429
mdlWorkerOut\(Att\)lmeInterceptCore 2 4 2843 2858 -1417 1 vs 2 22.785 < .001 
# Save variances
mdlWorkerOut$Att$varInterceptCore <-
  lme4::VarCorr(mdlWorkerOut$Att$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwInterceptCore <- 
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))
mdlWorkerOut$Att$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinInterceptCore <-
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))
mdlWorkerOut$Att$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))) * 100

We find that the the model with the added predictor indeed explains more variance in outgroup attitudes than the unconditional model. Looking at the individual coefficients, we find that the situational core need relates significantly to outgroup attitudes. The core need has little explained variance between participants (compared to the null model: Variance Explained = 1 – (Var with Predictor/Var without Predictor); -1.73%). The variance explained within participants is small to medium (7.21%).

random slope

In a next step, we check whether further letting the effect of core need fulfillment vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{9} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Att$lmeSlopesCore <-
  lme(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1,
    random = ~ 1 + keymotive_fulfillemt_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Att$lmerSlopesCore <- lmer(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 +
      (1 + keymotive_fulfillemt_1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2815.398
BIC 2839.148
Pseudo-R² (fixed effects) 0.046
Pseudo-R² (total) 0.763
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 71.863 66.133 77.594 24.578 19.326 0.000
keymotive_fulfillemt_1 0.179 0.061 0.297 2.970 18.186 0.008
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.111
PID keymotive_fulfillemt_1 0.219
Residual 7.941
Grouping Variables
Group # groups ICC
PID 21 0.732 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Att$lmeInterceptCore, 
      mdlWorkerOut$Att$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 5: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeInterceptCore 1 4 2843 2858 -1417
mdlWorkerOut\(Att\)lmeSlopesCore 2 6 2815 2839 -1402 1 vs 2 31.277 < .001 
# Save variances
mdlWorkerOut$Att$varSlopesCore <- 
  lme4::VarCorr(mdlWorkerOut$Att$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))
mdlWorkerOut$Att$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))
mdlWorkerOut$Att$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))) * 100

# Assumption Checks:
mdlWorkerOut$Att$diagSlopesCore <- 
  sjPlot::plot_model(mdlWorkerOut$Att$lmerSlopesCore, type = "diag")
grid.arrange(
  mdlWorkerOut$Att$diagSlopesCore[[1]],
  mdlWorkerOut$Att$diagSlopesCore[[2]]$`PID`,
  mdlWorkerOut$Att$diagSlopesCore[[3]],
  mdlWorkerOut$Att$diagSlopesCore[[4]]
)

# Plot prediction model
mdlWorkerOut$Att$predSlopesCore <- 
  workerOutgroupInteraction %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorkerOut$Att$lmeSlopesCore,
                           workerOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlWorkerOut$Att$predPltSlopesCore <-
    ggplot(data = mdlWorkerOut$Att$predSlopesCore, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/WorkerOut_PredictionPlot_SlopesAttCore.png",
  mdlWorkerOut$Att$predPltSlopesCore,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also find that the core need remains a strong predictor (even when letting the influence vary between participants).

TL;DR: The random slope adds significantly to the prediction model.

Need fulfillment and Interaction Quality

Based on the assertion that the relationship between core need fulfillment and outgroup attitudes is explained by a higher perceived interaction, the core need fulfillement should also significantly predict perceived interaction quality.

Unconditional model

Given that we now have the perceived interaction quality as our outcome variable of interest we again begin with an unconditional model (i.e., empty random intercept model), to see whether there is enough variance to explain within the participants. Similarly to before this is again done within the subsample of measurements during which an outgroup interaction was reported.

\[\begin{equation} \tag{10} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Qlt$lmerNull <-
  lme4::lmer(quality_overall_1 ~ 1 + (1 | PID), 
             data = workerOutgroupInteraction) # use optim if it does not converge
mdlWorkerOut$Qlt$lmeNull <-
  lme(
    quality_overall_1 ~ 1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.Qual.ML.r) #or with the lme function
summ(mdlWorkerOut$Qlt$lmerNull, digits = 3, center = TRUE)
Observations 387
Dependent variable quality_overall_1
Type Mixed effects linear regression
AIC 3347.534
BIC 3359.410
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.183
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 24.285 2.090 11.619 20.156 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 8.286
Residual 17.511
Grouping Variables
Group # groups ICC
PID 21 0.183 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1~1 + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlWorkerOut$Qlt$varNull <- 
  VarCorr(mdlWorkerOut$Qlt$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlWorkerOut$Qlt$tauNull <- 
  as.numeric(mdlWorkerOut$Qlt$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlWorkerOut$Qlt$sigmaNull <- 
  as.numeric(mdlWorkerOut$Qlt$varNull[2])
# The ICC estimate (between/between+within) is:
mdlWorkerOut$Qlt$IccNull <-
  (as.numeric(mdlWorkerOut$Qlt$varNull[1]) / (as.numeric(mdlWorkerOut$Qlt$varNull[1]) + as.numeric(mdlWorkerOut$Qlt$varNull[2])))
mdlWorkerOut$Qlt$IccPercNull <-
  ((as.numeric(mdlWorkerOut$Qlt$varNull[1]) / (as.numeric(mdlWorkerOut$Qlt$varNull[1]) + as.numeric(mdlWorkerOut$Qlt$varNull[2])))) * 100

We again find a reasonable amount of variance within the participants.

Note to self: For completeness, we should probably check this against an empty model without random intercept.

random intercept with level one predictor

We again add the core interaction need fulfillment to the multilevel random intercept model.

\[\begin{equation} \tag{11} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Qlt$lmeInterceptCore <-
  lme(
    quality_overall_1 ~ keymotive_fulfillemt_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(quality_overall_1 ~ keymotive_fulfillemt_1 + (1 | PID), 
       data = workerOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable quality_overall_1
Type Mixed effects linear regression
AIC 3287.252
BIC 3303.085
Pseudo-R² (fixed effects) 0.179
Pseudo-R² (total) 0.303
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 24.551 21.086 28.016 13.886 19.850 0.000
keymotive_fulfillemt_1 0.417 0.321 0.513 8.528 360.783 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 6.814
Residual 16.157
Grouping Variables
Group # groups ICC
PID 21 0.151 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(quality_overall_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Qlt$lmeNull, 
      mdlWorkerOut$Qlt$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 6: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Qlt\)lmeNull 1 3 3348 3359 -1671
mdlWorkerOut\(Qlt\)lmeInterceptCore 2 4 3287 3303 -1640 1 vs 2 62.283 < .001 
# Save variances
mdlWorkerOut$Qlt$varInterceptCore <-
  lme4::VarCorr(mdlWorkerOut$Qlt$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Qlt$varBtwInterceptCore <- 
  1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Qlt$varNull[1]))
mdlWorkerOut$Qlt$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Qlt$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Qlt$varWithinInterceptCore <-
  1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Qlt$varNull[2]))
mdlWorkerOut$Qlt$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Qlt$varNull[2]))) * 100

The predictor again adds a significant amount of explained variances beyond the empty model and looking at the slope coefficient, we find that the situational core need fulifillment relates significantly to perceived interaction quality. The core need explained substantial variance between participants (32.35%). The variance explained within participants is also medium (14.87%).

random slope

As before, we check whether further letting the effect of core need fulfillment vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{12} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Qlt$lmeSlopesCore <-
  lme(
    quality_overall_1 ~
      keymotive_fulfillemt_1,
    random = ~ 1 + keymotive_fulfillemt_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Qlt$lmerSlopesCore <-
    lmer(
      quality_overall_1 ~
        keymotive_fulfillemt_1 +
        (1 + keymotive_fulfillemt_1 | PID),
      data = workerOutgroupInteraction
    ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable quality_overall_1
Type Mixed effects linear regression
AIC 3289.190
BIC 3312.941
Pseudo-R² (fixed effects) 0.192
Pseudo-R² (total) 0.357
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 24.263 20.622 27.903 13.063 18.950 0.000
keymotive_fulfillemt_1 0.443 0.304 0.582 6.250 7.336 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 7.176
PID keymotive_fulfillemt_1 0.189
Residual 15.918
Grouping Variables
Group # groups ICC
PID 21 0.169 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(mdlWorkerOut$Qlt$lmerSlopesCore,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Qlt$lmeInterceptCore, 
      mdlWorkerOut$Qlt$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 7: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Qlt\)lmeInterceptCore 1 4 3287 3303 -1640
mdlWorkerOut\(Qlt\)lmeSlopesCore 2 6 3289 3313 -1639 1 vs 2 2.062 0.357 
# Save variances
mdlWorkerOut$Qlt$varSlopesCore <- 
  lme4::VarCorr(mdlWorkerOut$Qlt$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Qlt$varBtwSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]))
mdlWorkerOut$Qlt$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Qlt$varWithinSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]))
mdlWorkerOut$Qlt$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]))) * 100

We find that adding the random slopes does not add significantly beyond the random intercept model. This is unusual because this might indicate the the effect is very consistent across participants. However, we also see that when taking the possibility to varying slopes into account, the core need fulfillment remains a significant predictor of perceived interaction quality.

Need fulfillment, Interaction Quality, and Attitudes

In our final main step, we will jointly consider the effect of core need fulfillment and perceived interaction quality on outgroup attitudes. We expect that if the relationship between core need fulfillment and outgroup attitudes is indeed explained by a higher perceived interaction quality, the interaction quality perception should assume the explained variance of the core contact need fulfillment.

random intercept with level one predictors

We thus add both the core need fulfillment and perceived interaction quality to a random intercept model of outgroup attitudes.

\[\begin{equation} \tag{13} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmeInterceptCoreQlt <-
  lme(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + quality_overall_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + quality_overall_1 + (1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2756.572
BIC 2776.364
Pseudo-R² (fixed effects) 0.123
Pseudo-R² (total) 0.751
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 72.697 67.464 77.931 27.225 19.335 0.000
keymotive_fulfillemt_1 0.042 -0.009 0.094 1.614 371.091 0.107
quality_overall_1 0.252 0.204 0.300 10.288 367.215 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.025
Residual 7.574
Grouping Variables
Group # groups ICC
PID 21 0.716 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(
  mdlWorkerOut$Att$lmeNull, 
  mdlWorkerOut$Att$lmeInterceptCoreQlt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 8: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeNull 1 3 2863 2875 -1429
mdlWorkerOut\(Att\)lmeInterceptCoreQlt 2 5 2757 2776 -1373 1 vs 2 110.888 < .001 
# Save variances
mdlWorkerOut$Att$varInterceptCoreQlt <-
  lme4::VarCorr(mdlWorkerOut$Att$lmeInterceptCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwInterceptCoreQlt <- 
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))
mdlWorkerOut$Att$varBtwPercInterceptCoreQlt <- 
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinInterceptCoreQlt <-
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))
mdlWorkerOut$Att$varWithinPercInterceptCoreQlt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))) * 100

Unsurprisingly, the model with the two predictors adds significantly beyond the empty unconditional model. However, more importantly, looking at the coefficients, we find that the effect of core need fulfillemnt indeed is indeed strongly reduced and the variance is explained by the perceived interaction quality. The variance explained in outgroup attitudes is of medium effect size (between: 15.34%, within: 27.29%).

random slope

We again check whether further letting the effects vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{14} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

Check whether the omega matrix is correct should probably have as many taus as gammas.

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Att$lmeSlopesCoreQlt <-
  lme(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + quality_overall_1,
    random = ~ 1 + keymotive_fulfillemt_1 + quality_overall_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Att$lmerSlopesCoreQlt <- lmer(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + quality_overall_1 +
      (1 + keymotive_fulfillemt_1 + quality_overall_1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2669.434
BIC 2709.018
Pseudo-R² (fixed effects) 0.105
Pseudo-R² (total) 0.836
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 72.553 67.394 77.711 27.566 19.078 0.000
keymotive_fulfillemt_1 0.036 -0.076 0.149 0.636 15.027 0.535
quality_overall_1 0.233 0.128 0.338 4.345 21.118 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 11.827
PID keymotive_fulfillemt_1 0.216
PID quality_overall_1 0.210
Residual 6.133
Grouping Variables
Group # groups ICC
PID 21 0.788 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(
  mdlWorkerOut$Att$lmeNull,
  mdlWorkerOut$Att$lmeInterceptCoreQlt,
  mdlWorkerOut$Att$lmeSlopesCoreQlt
) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 9: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeNull 1 3 2863 2875 -1429
mdlWorkerOut\(Att\)lmeInterceptCoreQlt 2 5 2757 2776 -1373 1 vs 2 110.888 < .001
mdlWorkerOut\(Att\)lmeSlopesCoreQlt 3 10 2669 2709 -1325 2 vs 3 97.138 < .001 
# Save variances
mdlWorkerOut$Att$varSlopesCoreQlt <- 
  lme4::VarCorr(mdlWorkerOut$Att$lmeSlopesCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwSlopesCoreQlt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]))
mdlWorkerOut$Att$varBtwPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinSlopesCoreQlt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]))
mdlWorkerOut$Att$varWithinPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]))) * 100

# Assumption Checks:
mdlWorkerOut$Att$diagSlopesCoreQlt <- 
  sjPlot::plot_model(mdlWorkerOut$Att$lmerSlopesCoreQlt, type = "diag")
grid.arrange(
  mdlWorkerOut$Att$diagSlopesCoreQlt[[1]],
  mdlWorkerOut$Att$diagSlopesCoreQlt[[2]]$`PID`,
  mdlWorkerOut$Att$diagSlopesCoreQlt[[3]],
  mdlWorkerOut$Att$diagSlopesCoreQlt[[4]]
)

# Plot prediction model
mdlWorkerOut$Att$predSlopesCoreQlt <- 
  workerOutgroupInteraction %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorkerOut$Att$lmeSlopesCoreQlt,
                           workerOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlWorkerOut$Att$predPltSlopesCoreQlt <-
    ggplot(data = mdlWorkerOut$Att$predSlopesCoreQlt, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/WorkerOut_PredictionPlot_SlopesAttCoreQlt.png",
  mdlWorkerOut$Att$predPltSlopesCoreQlt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also find that the perceived interaction quality remains a strong predictor (even when letting the slopes vary between participants).

Check for robustness

To build further confidence in our results, we assess a few additional models that might offer alternative explanations of the effects we find.

Interaction Type

To make certain that the effect of core need fulfillment is specific to the interaction we compare the the effect to fulfillment of the situation core need when no outgroup interaction took place.

random intercept

Here we go back to the full dataset and add generalized situational core need fulfillment (either during an interaction or about the daytime in general) and whether an outgroup interaction happened as well as their interaction term to a random intercept model of outgroup attitudes.

\[\begin{equation} \tag{15} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorker$lmeInterceptAttCoreInt <-
  lme(
    thermometerDutch_1 ~ keyMotiveFulfilled * OutgroupInteraction,
    random =  ~ 1 | PID,
    data = workerInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlWorker$lmerInterceptAttCoreInt <- lmer(
    thermometerDutch_1 ~ keyMotiveFulfilled * OutgroupInteraction + (1 | PID),
    data = workerInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8774.952
BIC 8805.616
Pseudo-R² (fixed effects) 0.015
Pseudo-R² (total) 0.713
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.587 65.258 75.916 25.961 22.187 0.000
keyMotiveFulfilled 0.016 -0.007 0.039 1.377 1205.341 0.169
OutgroupInteraction 1.711 0.578 2.843 2.961 1202.807 0.003
keyMotiveFulfilled:OutgroupInteraction 0.109 0.061 0.157 4.428 1200.535 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.891
Residual 8.265
Grouping Variables
Group # groups ICC
PID 23 0.709 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeAttNull, 
      mdlWorker$lmeInterceptAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = c("l", rep("c", ncol(.)-1)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 10: Worker: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeAttNull 1 3 8806 8821 -4400
lmeInterceptAttCoreInt 2 6 8775 8806 -4381 1 vs 2 36.929 < .001 
# Save variances
mdlWorker$varInterceptAttCoreInt <- 
  lme4::VarCorr(mdlWorker$lmeInterceptAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwInterceptAttCoreInt <-
  1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[1]) / as.numeric(mdlWorker$varAttNull[1]))
mdlWorker$varBtwPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[1]) / as.numeric(mdlWorker$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinInterceptAttCoreInt <-
  1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[2]) / as.numeric(mdlWorker$varAttNull[2]))
mdlWorker$varWithinPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[2]) / as.numeric(mdlWorker$varAttNull[2]))) * 100

We find that the model explains significantly more variance than the empty null model. However, more interestingly, looking at the coefficients, we find that, as seen earlier, having an outgroup interaction has a strong effect on outgroup attitudes. Importantly, we find that there is no main effect of key need fulfillment but a significant interaction effect of core need fulfillment and outgroup contact. This indicates that it is not key need fulfillment in general — but only key need fulfillment during an outgroup contact that predicts more positive outgroup attitudes.

random slope

We again check whether further letting the effects vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{16} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorker$lmeSlopesAttCoreInt <- lme(
  thermometerDutch_1 ~
    keyMotiveFulfilled * OutgroupInteraction,
  random = ~ 1 + keyMotiveFulfilled + OutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = workerInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorker$lmerSlopesAttCoreInt <- lmer(
    thermometerDutch_1 ~
      keyMotiveFulfilled * OutgroupInteraction +
      (1 + keyMotiveFulfilled + OutgroupInteraction | PID),
    data = workerInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8763.482
BIC 8819.700
Pseudo-R² (fixed effects) 0.020
Pseudo-R² (total) 0.729
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.252 65.044 75.459 26.440 22.117 0.000
keyMotiveFulfilled 0.031 -0.021 0.082 1.166 10.093 0.270
OutgroupInteraction 1.846 0.169 3.523 2.157 19.514 0.044
keyMotiveFulfilled:OutgroupInteraction 0.112 0.057 0.167 3.984 396.330 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.558
PID keyMotiveFulfilled 0.100
PID OutgroupInteraction 2.795
Residual 8.043
Grouping Variables
Group # groups ICC
PID 23 0.709 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeAttNull, 
      mdlWorker$lmeInterceptAttCoreInt,
      mdlWorker$lmeSlopesAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 11: Worker: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeAttNull 1 3 8806 8821 -4400
lmeInterceptAttCoreInt 2 6 8775 8806 -4381 1 vs 2 36.929 < .001
lmeSlopesAttCoreInt 3 11 8763 8820 -4371 2 vs 3 21.47 < .001 
# Save variances
mdlWorker$varSlopesAttCoreInt <- 
  lme4::VarCorr(mdlWorker$lmeSlopesAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwSlopesAttCoreInt <- 
  1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[1]) / as.numeric(mdlWorker$varInterceptAttCoreInt[1]))
mdlWorker$varBtwPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[1]) / as.numeric(mdlWorker$varInterceptAttCoreInt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinSlopesAttCoreInt <- 
  1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[2]) / as.numeric(mdlWorker$varInterceptAttCoreInt[2]))
mdlWorker$varWithinPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[2]) / as.numeric(mdlWorker$varInterceptAttCoreInt[2]))) * 100

# Assumption Checks:
mdlWorker$diagSlopesAttCoreInt <-
  sjPlot::plot_model(mdlWorker$lmerSlopesAttCoreInt, type = "diag")
grid.arrange(
  mdlWorker$diagSlopesAttCoreInt[[1]],
  mdlWorker$diagSlopesAttCoreInt[[2]]$`PID`,
  mdlWorker$diagSlopesAttCoreInt[[3]],
  mdlWorker$diagSlopesAttCoreInt[[4]]
)

# Plot prediction model
mdlWorker$predSlopesAttCoreInt <- 
  workerInteractionType %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorker$lmeSlopesAttCoreInt,
                           workerInteractionType,
                           re.form = NA
                           )
         )

(
  mdlWorker$predPltSlopesAttCoreInt <-
    ggplot(data = mdlWorker$predSlopesAttCoreInt, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Worker_PredictionPlot_SlopesAttCoreInt.png",
  mdlWorker$predPltSlopesAttCoreInt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also see that when taking the possibility to varying slopes into account, the coefficient interpretations remains consistent (i.e., outgroup contact and its interaction with core need fulfillment remain significant predictors of positive outgroup attitudes).

Plot Interaction

Before we move on, we shortly illustrate the interaction effect of how the effect of core need fulfillment differed by whether an outgroup contact took place or not. To this end we illustrate (1) the raw data points (without taking the nested nature into account), as well as a plot of the model predicted values and their prediction interval (taking the nested structure into account based; similar to an interaction plot).

# visualize interaction
## Without ML structure
ggplot(data = workerInteractionType,
       aes(x = keyMotiveFulfilled,
           y = thermometerDutch_1,
           fill = OutgroupInteraction)) +
  #geom_point()+
  geom_smooth(method = 'lm',
              aes(linetype = OutgroupInteraction),
              color = "black") +
  #facet_wrap(~PID, ncol = 6)+
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  #scale_colour_manual(values=c("grey20", "black"), name="Intergroup Contact")+
  scale_y_continuous(
    limits = c(50, 100),
    breaks = seq(50, 100, by = 10),
    position = "left"
  ) +
  scale_x_continuous(limits = c(-50, 50), breaks = seq(-50, 50, by = 10)) +
  labs(
    title = "Without ML stucture",
    x = "Fulfillment Core Need",
    y = "Outgroup Attitudes",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact"
  ) +
  theme_Publication() +
  theme(legend.position = "bottom",
        legend.key.size = unit(1, "cm"))

## With ML structure
# create parameters for prediction
datNew = data.frame(
  OutgroupInteraction = factor(rep(c("No", "Yes"), each = 21)),
  keyMotiveFulfilled = rep(seq(-50, 50, 5), 2),
  PID = 0
)

# Predict values, clean up and calculate SE
PI <-
  merTools::predictInterval(
    merMod = mdlWorker$lmerSlopesAttCoreInt,
    newdata = datNew,
    level = 0.95,
    stat = "mean",
    type = "linear.prediction",
    include.resid.var = F,
    fix.intercept.variance = F
  )
mdlWorker$predInterceptAttCoreIntX <- 
  cbind(datNew, PI)
mdlWorker$predInterceptAttCoreIntX$se <-
  (mdlWorker$predInterceptAttCoreIntX$upr - mdlWorker$predInterceptAttCoreIntX$fit) / 1.96
rm(datNew, PI)

# Plot predicted values with SE
ggplot(
  mdlWorker$predInterceptAttCoreIntX,
  aes(x = keyMotiveFulfilled,
      y = fit,
      fill = OutgroupInteraction)
)+
  #geom_point() +
  geom_line(aes(linetype = as.factor(OutgroupInteraction)), size = 1) +
  #facet_wrap(~PID, ncol = 6)+
  geom_ribbon(data = mdlWorker$predInterceptAttCoreIntX,
              aes(ymin = fit - se, ymax = fit + se),
              alpha = 0.3) +
  scale_x_continuous(breaks = seq(-50, 50, 10)) +
  scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  labs(
    x = "Fulfillment Core Need",
    y = "Outgroup Attitude (NL)",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact",
    title = "Based on Model Predictions"
  ) +
  theme_Publication()

# #### Bayesian estimation !! ONLY RUN ON FINAL RENDER !! Takes forever ####
# options(mc.cores = parallel::detectCores())  # Run many chains simultaneously
# brmfit <- brm(
#   thermometerDutch_1 ~ keyMotiveFulfilled * OutgroupInteraction + 
#     (1 + keyMotiveFulfilled + OutgroupInteraction | PID),
#   data = workerInteractionType,
#   family = gaussian,
#   iter = 1000,
#   chains = 4
# )
# 
# # create parameters for prediction:
# datNew = data.frame(
#   OutgroupInteraction = rep(c("No", "Yes"), each = 101),
#   keyMotiveFulfilled = rep(seq(-50, 50, 1), 2)
# )
# # Save predicted values and adjust names and labels
# fitavg <-
#   cbind(datNew,
#         fitted(brmfit, newdata = datNew, re_formula = NA)[, -2])
# names(fitavg)[names(fitavg) == "Estimate"] = "pred"
# fitavg$se <- (fitavg$Q97.5 - fitavg$pred) / 1.96
# 
# # Plot Bayesian SE prediction interval
# ggplot(fitavg,
#        aes(x = keyMotiveFulfilled,
#            y = pred,
#            fill = OutgroupInteraction)) +
#   scale_x_continuous(breaks = seq(-50, 50, 10)) +
#   scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
#   geom_line(aes(linetype = OutgroupInteraction), size = 1) +
#   geom_ribbon(aes(ymin = pred - se, ymax = pred + se), alpha = 0.2) +
#   scale_linetype_manual(values = c("dashed", "solid")) +
#   scale_fill_manual(values = c("darkgrey", "black")) +
#   labs(
#     x = "Fulfillment Core Need",
#     y = "Outgroup Attitude (NL)",
#     fill = "Intergroup Contact",
#     linetype = "Intergroup Contact",
#     title = "Based on Bayesian Prediction Interval"
#   ) +
#   theme_Publication()
# 
# # plot all overlayed posteriors:
# pst <- posterior_samples(brmfit, "b")
# ggplot(workerInteractionType,
#        aes(x = keyMotiveFulfilled, y = thermometerDutch_1)) +
#   geom_point(shape = 4, alpha = .1) +
#   geom_tile() +
#   geom_abline(
#     data = pst,
#     aes(intercept = b_Intercept, slope = b_keyMotiveFulfilled),
#     alpha = .025,
#     size = .4
#   ) +
#   labs(title = "slope Posteriors",
#        x = "Fulfillment Core Need",
#        y = "Outgroup Attitudes (NL)") +
#   theme_Publication()
# rm(datNew, brmfit, fitavg, pst)

Other psychological needs

In a final step we check whether during the interaction the core situational need is a meaningful predictor even when taking other fundamental psychological needs into account. We focus on the three commonly considered self determination needs: competence, autonomy, and relatedness.

random intercept with level ome predictors

We add the core need fulfillment with the three self determination needs to a random intercept model of outgroup attitudes.

\[\begin{equation} \tag{17} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmeInterceptCoreSdt <-
  lme(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction,
    na.action = na.exclude
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlWorkerOut$Att$lmerInterceptCoreSdt <- lmer(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 + (1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 386
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2804.147
BIC 2831.838
Pseudo-R² (fixed effects) 0.086
Pseudo-R² (total) 0.754
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 73.077 67.364 78.790 25.071 19.379 0.000
keymotive_fulfillemt_1 0.092 0.038 0.146 3.323 365.656 0.001
competence_1 0.019 -0.034 0.072 0.705 363.428 0.481
autonomy_1 0.079 0.014 0.143 2.388 363.976 0.017
relatedness_1 0.114 0.078 0.150 6.204 365.509 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.139
Residual 7.978
Grouping Variables
Group # groups ICC
PID 21 0.731 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# To be compared against a model with only the self determination theory needs
mdlWorkerOut$Att$lmeInterceptSdt <-
  lme(
    thermometerDutch_1 ~ competence_1 + autonomy_1 + relatedness_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction,
    na.action = na.exclude
  )

summ(
  mdlWorkerOut$Att$lmerInterceptSdt <- lmer(
    thermometerDutch_1 ~ competence_1 + autonomy_1 + relatedness_1 + (1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 386
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2807.710
BIC 2831.445
Pseudo-R² (fixed effects) 0.073
Pseudo-R² (total) 0.745
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 73.036 67.331 78.741 25.091 19.368 0.000
competence_1 0.041 -0.011 0.093 1.546 365.717 0.123
autonomy_1 0.092 0.027 0.157 2.778 365.087 0.006
relatedness_1 0.119 0.082 0.155 6.382 366.680 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.114
Residual 8.089
Grouping Variables
Group # groups ICC
PID 21 0.724 
# Compare new model to previous steps
anova(
  mdlWorkerOut$Att$lmeInterceptSdt,
  mdlWorkerOut$Att$lmeInterceptCoreSdt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 12: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeInterceptSdt 1 6 2808 2831 -1398
mdlWorkerOut\(Att\)lmeInterceptCoreSdt 2 7 2804 2832 -1395 1 vs 2 5.563 0.018 
rm(lmeInterceptCoreRed)

# Save variances
mdlWorkerOut$Att$varInterceptCoreSdt <-
  lme4::VarCorr(mdlWorkerOut$Att$lmeInterceptCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwInterceptCoreSdt <- 
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))
mdlWorkerOut$Att$varBtwPercInterceptCoreSdt <- 
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinInterceptCoreSdt <-
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))
mdlWorkerOut$Att$varWithinPercInterceptCoreSdt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))) * 100

We compare the models of the core need and the SDT need fulfillments to a model that only includes the SDT needs. We find that the core need adds significantly above the SDT needs. We find that next to relatedness, the core need explains the most variance and compared to the model with only the SDT needs, the core need fulfillment flexibly takes on some of the explained variance of all of the three fundamental needs (i.e., reduction in SDT beta weights).

random slope

We again check whether further letting the effects vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{18} Attitude_{ti} = \\gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Att$lmeSlopesCoreSdt <-
  lme(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1,
    random = ~ 1 + keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction,
    na.action = na.exclude
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Att$lmerSlopesCoreSdt <- lmer(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 +
      (1 + keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 386
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2754.434
BIC 2837.507
Pseudo-R² (fixed effects) 0.067
Pseudo-R² (total) 0.870
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.956 65.168 76.745 24.027 18.872 0.000
keymotive_fulfillemt_1 0.091 0.012 0.169 2.259 10.955 0.045
competence_1 0.047 -0.166 0.259 0.429 9.300 0.678
autonomy_1 0.059 -0.120 0.239 0.646 9.357 0.534
relatedness_1 0.100 0.056 0.145 4.412 18.715 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.289
PID keymotive_fulfillemt_1 0.123
PID competence_1 0.467
PID autonomy_1 0.378
PID relatedness_1 0.064
Residual 6.423
Grouping Variables
Group # groups ICC
PID 21 0.811 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Att$lmeInterceptCoreSdt, 
      mdlWorkerOut$Att$lmeSlopesCoreSdt) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 13: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeInterceptCoreSdt 1 7 2804 2832 -1395
mdlWorkerOut\(Att\)lmeSlopesCoreSdt 2 21 2756 2838 -1357 1 vs 2 76.506 < .001 
# Save variances
mdlWorkerOut$Att$varSlopesCoreSdt <- 
  lme4::VarCorr(mdlWorkerOut$Att$lmeSlopesCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwSlopesCoreSdt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]))
mdlWorkerOut$Att$varBtwPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinSlopesCoreSdt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]))
mdlWorkerOut$Att$varWithinPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]))) * 100

# Assumption Checks:
mdlWorkerOut$Att$diagSlopesCoreSdt <- 
  sjPlot::plot_model(mdlWorkerOut$Att$lmerSlopesCoreSdt, type = "diag")
grid.arrange(
  mdlWorkerOut$Att$diagSlopesCoreSdt[[1]],
  mdlWorkerOut$Att$diagSlopesCoreSdt[[2]]$`PID`,
  mdlWorkerOut$Att$diagSlopesCoreSdt[[3]],
  mdlWorkerOut$Att$diagSlopesCoreSdt[[4]]
)

# Plot prediction model
mdlWorkerOut$Att$predSlopesCoreSdt <- 
  workerOutgroupInteraction %>%
  filter(!is.na(autonomy_1)) %>%
  select(thermometerDutch_1, session, PID, autonomy_1) %>% 
  mutate(measure = predict(mdlWorkerOut$Att$lmeSlopesCoreSdt,
                           workerOutgroupInteraction %>% filter(!is.na(autonomy_1)),
                           re.form = NA
                           )
         )

(
  mdlWorkerOut$Att$predPltSlopesCoreSdt <-
    ggplot(data = mdlWorkerOut$Att$predSlopesCoreSdt, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/WorkerOut_PredictionPlot_SlopesAttCoreStd.png",
  mdlWorkerOut$Att$predPltSlopesCoreSdt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also see that when taking the possibility to varying slopes into account, the coefficient interpretations remains consistent (i.e., core need and relatedness remain the strongest and only significant predictors). Note however a slight reduction in the p value of the core need fulfillment.

Student Sample

The second sample we assess is larger study with international psychology students at Western European university. The hypotheses mirror those of the first study and are re-iterated here:

  1. Based on the most general understanding of the contact hypothesis, an increase in frequency and quality of contact should jointly account for changes in more favorable outgroup attitudes.
    1. Participants with more intergroup interactions should have a more favorable outgroup attitudes.
    2. Outgroup attitudes should be higher after an intergroup interaction compared to a non-outgroup interaction.
    3. Participants with more intergroup interactions should have a more favorable outgroup attitudes depending on the average interaction quality.
  2. Based on our proposal, intergroup interactions with higher situational core need fulfillment should predict more favorable outgroup attitudes due to more positive interaction quality perceptions.
    1. Outgroup attitudes should be more favorable after intergroup interactions with high key need fulfillment.
    2. Interaction Quality should be perceived as more positive after intergroup interactions with higher key need fulfillment.
    3. The variance explained in outgroup attitudes by key need fulfillment should to a large extend be assumed by interaction quality.
    4. The effect of key need fulfillment on outgroup attitudes should be specific to intergroup interactions and not be due to need fulfillment in general. Thus, the effect of key need fulfillment on outgroup attitudes should stronger for intergroup interact than for ingroup interactions.
    5. The effect of key need fulfillment on outgroup attitudes should be persist even when taking other fundamental psychological needs into account. Thus, the effect of key need fulfillment on outgroup attitudes should remain strong even after controlling for autonomy, competence, and relatedness fulfillment during the interaction (cf., self-determination theory).

We will test our main hypotheses for this study in a sequential manner.

Data Description

Still in ‘scr/studentDescriptive.R’. Needs to be merged with this document.

Contact Hypothesis

We again test the most general contact hypothesis in two steps. First, we assess whether more intergroup interactions are related to to more positive outgroup attitudes. Second, we test whether a potential positive effect on outgroup attitudes depends on the interaction quality (jointly with the number of interactions).

Interaction Frequency and Outgroup Attitudes

To test the impact of the overall number of interactions, we hope to find a significant relationship between the number of interactions a participant had and the average outgroup attitude.

\[\begin{equation} \tag{19} r_{ContactFrequency, OutgroupAttitudes} \neq 0 \end{equation}\]

# correlation panel
pairs.panels.new(
  studentContactFreq %>% select(SumContactNL, SumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  studentContactFreq %>% select(WinSumContactNL, WinSumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

We find that both the number of interactions and the number of measurement beeps with an interaction are significantly related with the average outgroup attitudes. This is to say that within our data, participants with more outgroup interactions did have significantly more positive outgroup attitudes. This is inconsistent with the results we found in the worker sample.

Outgroup Attitudes by Interaction Type

In a next step we take into account that having an interaction with an outgroup member, does not happen in a social vacuum. Participants who indicated that they had an interaction with an outgroup member include measurement occasions during which someone either only had an interaction with an outgroup member as well as times during which a person also had interaction(s) with a non-Dutch person. Inversely, participants who indicated that they did not have an interaction with a Dutch person might either have had no interaction at all or had an interaction with a non-Dutch person. We, thus consider all possible combinations and their independent influences on outgroup attitudes.

OLS regression

We first assess the impact of the different interaction types across all measurement points (lumping all beeps together).

\[\begin{equation} \tag{20} \mu_{i,0utgroupInteraction} > \mu_{i,IngroupInteraction} \\ Attitude = OutgroupInteraction \times NonOutgroupInteraction \end{equation}\]

# between participants interaction type
studentContactType <- studentInteractionType %>%
  group_by(
    OutgroupInteraction,
    NonOutgroupInteraction
  ) %>%
  summarise(
    n = n(),
    AttitudeM = mean(AttitudesDutch, na.rm = TRUE),
    AttitudeSD = sd(AttitudesDutch, na.rm = TRUE),
    AttitudeSE = AttitudeSD / sqrt(n),
    AttitudeLwr = AttitudeM - 1.96 * AttitudeSE,
    AttitudeUpr = AttitudeM + 1.96 * AttitudeSE
  ) %>%
  ungroup()

# plot bar chart
ggplot(
  studentContactType,
  aes(
    y = AttitudeM,
    x = OutgroupInteraction,
    fill = NonOutgroupInteraction
  )
) +
  geom_bar(
    stat = "identity",
    color = "black",
    position = position_dodge()
  ) +
  geom_errorbar(aes(ymin = AttitudeM, ymax = AttitudeUpr),
    width = .2,
    position = position_dodge(.9)
  ) +
  labs(
    fill = "Non-Outgroup Interaction",
    x = "Outgroup Interaction",
    y = "Outgroup Attitude",
    title = "Outgroup Attitudes by Interaction Type [95% CI]"
  ) +
  scale_fill_grey(
    start = 0.2,
    end = 0.8
  ) +
  scale_y_continuous(
    limits = c(0, 100),
    breaks = c(0, 15, 30, 40, 50, 60, 70, 85, 100),
    labels = c(
      "Very cold or unfavorable feelings 0°",
      "Quite cold and unfavorable feelings 15°",
      "Fairly cold and unfavorable feelings 30°",
      "A bit cold and unfavorable feelings 40°",
      "No feeling at all 50°",
      "A bit warm and favorable feelings 60°",
      "Fairly warm and favorable feelings 70° ",
      "Quite warm and favorable feelings 85° ",
      "Very warm and favorable feelings 100° "
    )
  ) +
  theme_Publication()

# create list to store student models
mdlStudent <- list()

# regression
mdlStudent$lmAttInt <-
  lm(AttitudesDutch ~ OutgroupInteraction * NonOutgroupInteraction,
    data = studentInteractionType
  )
# summary(lmstudentAttInteraction)

summ(
  mdlStudent$lmAttInt,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type OLS linear regression
F(3,4961) 37.435
R² 0.022
Adj. R² 0.022
Est. 2.5% 97.5% t val. p
(Intercept) 66.807 65.815 67.799 132.053 0.000
OutgroupInteraction 8.159 5.699 10.619 6.503 0.000
NonOutgroupInteraction -1.226 -2.529 0.077 -1.845 0.065
OutgroupInteraction:NonOutgroupInteraction -0.355 -3.439 2.730 -0.225 0.822
Standard errors: OLS; Continuous predictors are mean-centered.

We find that while controlling for interactions with non-Dutch people, outgroup attitudes were significantly higher when participants had an interaction with a Dutch person. The effect is of a medium size (8.16 points on a 0–100 scale). However, this analysis lumps all ESM beeps from every participants together and ignores that the data is nested within participants.

ML Regression with interaction term

We, thus, proceed with a multilevel analysis, which acknowledges that the measurements are nested within participants.

Unconditional model

We start by creating an empty random intercept model (i.e., let the outgroup attitude intercept be different between participants; unconditional model).

\[\begin{equation} \tag{21} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudent$lmerAttNull <-
  lme4::lmer(AttitudesDutch ~ 1 + (1 | PID),
    data = dtStudents$full
  ) # use optim if it does not converge

mdlStudent$lmeAttNull <-
  lme(
    AttitudesDutch ~ 1,
    random = ~ 1 | PID,
    data = dtStudents$full,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(mdlStudent$lmerAttNull) #or with the lme function
summ(mdlStudent$lmerAttNull, digits = 3, center = TRUE)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36740.926
BIC 36760.456
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.801
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 67.290 1.751 38.435 111.704 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.540
Residual 9.235
Grouping Variables
Group # groups ICC
PID 113 0.801 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch~1 + (1|PID),data=dtStudent$full),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlStudent$varAttNull <- 
  VarCorr(mdlStudent$lmeAttNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlStudent$tauAttNull <- 
  as.numeric(mdlStudent$varAttNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlStudent$sigmaAttNull <- 
  as.numeric(mdlStudent$varAttNull[2])
# The ICC estimate (between/between+within) is:
mdlStudent$IccAttNull <-
  (as.numeric(mdlStudent$varAttNull[1]) / (as.numeric(mdlStudent$varAttNull[1]) + as.numeric(mdlStudent$varAttNull[2])))
mdlStudent$IccPercAttNull <-
  ((as.numeric(mdlStudent$varAttNull[1]) / (as.numeric(mdlStudent$varAttNull[1]) + as.numeric(mdlStudent$varAttNull[2])))) * 100

We then compare the random intercept model to a model without a random intercept (i.e., without levels at all).

# Create and save Model
mdlStudent$glsAttNull  <-
  gls(AttitudesDutch ~ 1,
      data = dtStudents$full,
      control = list(opt = "nlmimb"))

# calculate Deviances manually:
mdlStudent$DevianceGlsNull <- logLik(mdlStudent$glsAttNull) * -2
mdlStudent$DevianceMlNull <- logLik(mdlStudent$lmeAttNull) * -2

# Compare the two null models:
anova(mdlStudent$glsAttNull,
      mdlStudent$lmeAttNull) %>% 
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 14: Student: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
glsAttNull 1 2 44352 44365 -22174
lmeAttNull 2 3 36741 36760 -18367 1 vs 2 7613.538 < .001

Comparing the two empty model, we find that there is indeed a significant amount of variance explained by including a random intercept.

To assess the variances within and between participants we look at the ICC \(\tau_{00}^2 / (\tau_{00}^2 + \sigma^2)\): The ratio of the between-cluster variance to the total variance is called the Intraclass Correlation. It tells you the proportion of the total variance in Y that is accounted for by the clustering. (In this case the clustering means clustering observations per participant).

We find that an estimated 80.12% of the variation in Feeling Thermometer scores is explained by between participant differences (clustering by PID). This is to say that 80.12% of the variance in any individual report of Attitudes towards the Dutch can be explained by the properties of the individual who provided the rating. And we find that including ‘participant’ as a predictor adds significantly to the model.

random intercept with predictors

To this random intercept model we now add the two types of interactions possible at each measurement point as contemporaneous predictors of outgroup attitudes. That is: We check whether within participants having an outgroup interaction (or a non-outgroup interaction) is associated with more positive outgroup attitudes at the same measurement point.

\[\begin{equation} \tag{22} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudent$lmeInterceptAttType <-
  lme(
    AttitudesDutch ~ OutgroupInteraction + NonOutgroupInteraction,
    random =  ~ 1 | PID,
    data = studentInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    AttitudesDutch ~ OutgroupInteraction + NonOutgroupInteraction + (1 | PID),
    data = studentInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36684.392
BIC 36716.943
Pseudo-R² (fixed effects) 0.003
Pseudo-R² (total) 0.801
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.800 63.374 70.226 38.217 114.319 0.000
OutgroupInteraction 2.882 2.162 3.601 7.852 4862.093 0.000
NonOutgroupInteraction -0.133 -0.708 0.443 -0.453 4861.600 0.651
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.400
Residual 9.180
Grouping Variables
Group # groups ICC
PID 113 0.801 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentInteractionType),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeAttNull, 
      mdlStudent$lmeInterceptAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 15: Student: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeAttNull 1 3 36741 36760 -18367
lmeInterceptAttType 2 5 36684 36717 -18337 1 vs 2 60.533 < .001 
# Save variances
mdlStudent$varInterceptAttType <- 
  lme4::VarCorr(mdlStudent$lmeInterceptAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwInterceptAttType <-
  1 - (as.numeric(mdlStudent$varInterceptAttType[1]) / as.numeric(mdlStudent$varAttNull[1]))
mdlStudent$varBtwPercInterceptAttType <-
  (1 - (as.numeric(mdlStudent$varInterceptAttType[1]) / as.numeric(mdlStudent$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinInterceptAttType <-
  1 - (as.numeric(mdlStudent$varInterceptAttType[2]) / as.numeric(mdlStudent$varAttNull[2]))
mdlStudent$varWithinPercInterceptAttType <-
  (1 - (as.numeric(mdlStudent$varInterceptAttType[2]) / as.numeric(mdlStudent$varAttNull[2]))) * 100

We find that a random intercept model with the two interaction types as predictors explains significantly more variance then an empty random intercept model. Looking at the individual coefficients, we find that having an outgroup interaction is indeed associated with significantly more positive outgroup attitudes, while having an interaction with a non-Dutch person does not significantly relate to more positive or negative attitudes towards the Dutch.

TL;DR: Interaction with Dutch is great predictor, interactions with non-Dutch is not.

random slope

In a next step, we check whether further letting the effect of the different interaction types vary between participants explains additional variance in outgroup attitudes (i.e., random slopes).

\[\begin{equation} \tag{5} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudent$lmeSlopesAttType <- lme(
  AttitudesDutch ~
    OutgroupInteraction + NonOutgroupInteraction,
  random = ~ 1 + OutgroupInteraction + NonOutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = studentInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudent$lmerSlopesAttType <- lmer(
    AttitudesDutch ~
      OutgroupInteraction + NonOutgroupInteraction +
      (1 + OutgroupInteraction + NonOutgroupInteraction | PID),
    data = studentInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36476.495
BIC 36541.597
Pseudo-R² (fixed effects) 0.003
Pseudo-R² (total) 0.817
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.717 63.293 70.141 38.191 111.810 0.000
OutgroupInteraction 2.992 1.448 4.537 3.797 93.820 0.000
NonOutgroupInteraction 0.009 -0.763 0.781 0.022 108.800 0.982
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.375
PID OutgroupInteraction 6.996
PID NonOutgroupInteraction 2.689
Residual 8.801
Grouping Variables
Group # groups ICC
PID 113 0.813 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeInterceptAttType, 
      mdlStudent$lmeSlopesAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 16: Student: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeInterceptAttType 1 5 36684 36717 -18337
lmeSlopesAttType 2 10 36476 36542 -18228 1 vs 2 217.896 < .001 
# Save variances
mdlStudent$varSlopesAttType <- 
  lme4::VarCorr(mdlStudent$lmeSlopesAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwSlopesAttType <- 
  1 - (as.numeric(mdlStudent$varSlopesAttType[1]) / as.numeric(mdlStudent$varInterceptAttType[1]))
mdlStudent$varBtwPercSlopesAttType <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttType[1]) / as.numeric(mdlStudent$varInterceptAttType[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinSlopesAttType <- 
  1 - (as.numeric(mdlStudent$varSlopesAttType[2]) / as.numeric(mdlStudent$varInterceptAttType[2]))
mdlStudent$varWithinPercSlopesAttType <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttType[2]) / as.numeric(mdlStudent$varInterceptAttType[2]))) * 100

# Assumption Checks:
mdlStudent$diagSlopesAttType <-
  sjPlot::plot_model(mdlStudent$lmerSlopesAttType, type = "diag")
grid.arrange(
  mdlStudent$diagSlopesAttType[[1]],
  mdlStudent$diagSlopesAttType[[2]]$`PID`,
  mdlStudent$diagSlopesAttType[[3]],
  mdlStudent$diagSlopesAttType[[4]]
)

# Plot prediction model
mdlStudent$predSlopesAttType <- 
  studentInteractionType %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudent$lmeSlopesAttType,
                           studentInteractionType,
                           re.form = NA
                           )
         )

(
  mdlStudent$predPltSlopesAttType <-
    ggplot(data = mdlStudent$predSlopesAttType %>% filter(PID %in% studentPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Student_PredictionPlot_SlopesAttType.png",
  mdlStudent$predPltSlopesAttType,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. This is is different to the previous study where the random slope did not add significantly.

TL;DR: Random slopes adds significantly.

Interaction Frequency and Interaction Quality

In a final step we check whether the effect outgroup interactions, in part, depends on the quality during the interactions. Because we can only assess interaction quality when there is an interaction, it is difficult to assess this with the interaction dummy as a within person predictor. Instead, we will use an aggregate measure of interaction quality and average interaction quality to consider the two predictors jointly.

\[\begin{equation} \tag{6} Attitude = ContactFreq \times AverageContactQual \end{equation}\]

# correlation panel
pairs.panels.new(
  studentContactFreq %>% select(SumContactNL, SumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  studentContactFreq %>% select(WinSumContactNL, WinSumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

Within the data, we find no significant correlation between the participants’ Average Interaction Quality and their Average Outgroup Attitudes. Thus, within our data participants with a higher quality outgroup interactions did not hold more positive attitudes towards that group. However, the frequency of intergroup interactions had a meaningful correlation with both the average interaction quality or average outgroup attitudes.

# regression
mdlStudent$lmAttFreqQualX <-
  lm(AvAttitude ~ SumContactNL * AvQuality, data = studentContactFreq)
# summary(mdlStudent$lmAttFreqQualX)

summ(
  mdlStudent$lmAttFreqQualX,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 113
Dependent variable AvAttitude
Type OLS linear regression
F(3,109) 3.748
R² 0.094
Adj. R² 0.069
Est. 2.5% 97.5% t val. p
(Intercept) 67.618 64.156 71.080 38.707 0.000
SumContactNL 0.805 0.320 1.290 3.291 0.001
AvQuality 0.253 -0.135 0.641 1.292 0.199
SumContactNL:AvQuality 0.023 -0.032 0.078 0.819 0.414
Standard errors: OLS; Continuous predictors are mean-centered. 
interactions::interact_plot(
  mdlStudent$lmAttFreqQualX,
  pred = SumContactNL,
  modx = AvQuality,
  interval = TRUE,
  partial.residuals = TRUE
)

interactions::johnson_neyman(mdlStudent$lmAttFreqQualX,
                             pred = SumContactNL,
                             modx = AvQuality,
                             alpha = .05)
## JOHNSON-NEYMAN INTERVAL 
## 
## When AvQuality is INSIDE the interval [73.80, 100.18], the slope of SumContactNL is p < .05.
## 
## Note: The range of observed values of AvQuality is [61.50, 100.00]

We find that in our student sample there is only a relationship between the number of outgroup contacts but no significant effect of average perceived contact quality. Nor do we find that in this sample the impact of the number of interactions is moderated by the average contact quality. This is not entirely consistent with the sojourner sample, where average contact quality did have a meaningful effect on outgroup attitudes. This effect is not necessarily surprising given that the variables aggregate all within person variation and there were substantially more measurements where participants did not have an interaction (but reported their outgroup attitudes) than measurements that followed an outgroup contact.

Need Fulfillment

The main proposal of our article is that the success of an outgroup interaction might be explained by whether or not the interaction fulfilled the person’s core situational need. This should, in turn, be due to a higher perceived interaction quality. We will this sequentially test whether the fulfillment of the core need during an interaction is (1) related to more positive outgroup attitudes, (2) higher perceived interaction quality, and (3) whether the variance explained by the core need is assumed by the perceived interaction quality if considered jointly.

Need fulfillment and Attitudes

In a first step we, thus, test the relationship between outgroup attitudes and the fulfillment of the core situational need during the interaction.

Unconditional model

We again start by creating an empty random intercept model (i.e., let the outgroup attitude intercept be different between participants; unconditional model). Note that this unconditional model differs from the empty model created earlier because for this set of analyses we will only consider the subsample of measurement points during which an outgroup interaction was reported. This is necessary because measurements of needs during the interaction and perceived interaction quality are only meaningful within an interaction context.

# see how large our outgroup interaction subset actually is
tbl_cross(
  studentInteractionType,
  row = OutgroupInteraction,
  col = NonOutgroupInteraction,
  percent = "row"
)
Characteristic NonOutgroupInteraction Total
No Yes
OutgroupInteraction
No 1,695 (42%) 2,335 (58%) 4,030 (100%)
Yes 329 (35%) 606 (65%) 935 (100%)
Total 2,024 (41%) 2,941 (59%) 4,965 (100%)
 
# create outgroup interaction subset
studentOutgroupInteraction <- studentInteractionType %>%
  filter(OutgroupInteraction == "Yes")

# create empty list to organize models
mdlStudentOut <- 
  list(
    Att = list(),
    Qlt = list()
  )

Note to self: For completeness, we should probably check this against an empty model without random intercept.

\[\begin{equation} \tag{23} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmerNull <-
  lme4::lmer(AttitudesDutch ~ 1 + (1 | PID), 
             data = studentOutgroupInteraction) # use optim if it does not converge
mdlStudentOut$Att$lmeNull <-
  lme(
    AttitudesDutch ~ 1,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.ML.r) #or with the lme function
summ(mdlStudentOut$Att$lmerNull, digits = 3, center = TRUE)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7250.070
BIC 7264.591
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.724
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 70.736 1.609 43.954 102.806 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 16.042
Residual 9.898
Grouping Variables
Group # groups ICC
PID 108 0.724 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch~1 + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlStudentOut$Att$varNull <- 
  VarCorr(mdlStudentOut$Att$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlStudentOut$Att$tauNull <- 
  as.numeric(mdlStudentOut$Att$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlStudentOut$Att$sigmaNull <- 
  as.numeric(mdlStudentOut$Att$varNull[2])
# The ICC estimate (between/between+within) is:
mdlStudentOut$Att$IccNull <-
  (as.numeric(mdlStudentOut$Att$varNull[1]) / (as.numeric(mdlStudentOut$Att$varNull[1]) + as.numeric(mdlStudentOut$Att$varNull[2])))
mdlStudentOut$Att$IccPercNull <-
  ((as.numeric(mdlStudentOut$Att$varNull[1]) / (as.numeric(mdlStudentOut$Att$varNull[1]) + as.numeric(mdlStudentOut$Att$varNull[2])))) * 100

random intercept with level one predictors

We thus add the core interaction need fulfillment to the multilevel random intercept model.

\[\begin{equation} \tag{8} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmeInterceptCore <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(AttitudesDutch ~ KeyNeedFullfillment + (1 | PID), 
       data = studentOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7234.687
BIC 7254.049
Pseudo-R² (fixed effects) 0.009
Pseudo-R² (total) 0.728
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.811 67.685 73.936 44.404 102.944 0.000
KeyNeedFullfillment 0.096 0.057 0.134 4.866 851.477 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.900
Residual 9.777
Grouping Variables
Group # groups ICC
PID 108 0.726 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Att$lmeNull, 
      mdlStudentOut$Att$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 17: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeNull 1 3 7250 7265 -3622
mdlStudentOut\(Att\)lmeInterceptCore 2 4 7235 7254 -3613 1 vs 2 17.383 < .001 
# Save variances
mdlStudentOut$Att$varInterceptCore <-
  lme4::VarCorr(mdlStudentOut$Att$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwInterceptCore <- 
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))
mdlStudentOut$Att$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinInterceptCore <-
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))
mdlStudentOut$Att$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))) * 100

We find that the the model with the added predictor indeed explains more variance in outgroup attitudes than the unconditional model. Looking at the individual coefficients, we find that the situational core need relates significantly to outgroup attitudes. The core need has little explained variance between participants (compared to the null model: Variance Explained = 1 – (Var with Predictor/Var without Predictor); 1.77%). The variance explained within participants is small to medium (2.44%).

random slope

In a next step, we check whether further letting the effect of core need fulfillment vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{9} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudentOut$Att$lmeSlopesCore <-
  lme(
    AttitudesDutch ~
      KeyNeedFullfillment,
    random = ~ 1 + KeyNeedFullfillment | PID,
    control = lmeControl(opt = "optim"),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Att$lmerSlopesCore <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment +
      (1 + KeyNeedFullfillment | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7218.921
BIC 7247.965
Pseudo-R² (fixed effects) 0.017
Pseudo-R² (total) 0.745
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.915 67.809 74.022 44.739 102.943 0.000
KeyNeedFullfillment 0.132 0.073 0.191 4.382 45.436 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.774
PID KeyNeedFullfillment 0.150
Residual 9.472
Grouping Variables
Group # groups ICC
PID 108 0.735 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Att$lmeInterceptCore, 
      mdlStudentOut$Att$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 18: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeInterceptCore 1 4 7235 7254 -3613
mdlStudentOut\(Att\)lmeSlopesCore 2 6 7219 7248 -3603 1 vs 2 19.759 < .001 
# Save variances
mdlStudentOut$Att$varSlopesCore <- 
  lme4::VarCorr(mdlStudentOut$Att$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))
mdlStudentOut$Att$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))
mdlStudentOut$Att$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))) * 100

# Assumption Checks:
mdlStudentOut$Att$diagSlopesCore <- 
  sjPlot::plot_model(mdlStudentOut$Att$lmerSlopesCore, type = "diag")
grid.arrange(
  mdlStudentOut$Att$diagSlopesCore[[1]],
  mdlStudentOut$Att$diagSlopesCore[[2]]$`PID`,
  mdlStudentOut$Att$diagSlopesCore[[3]],
  mdlStudentOut$Att$diagSlopesCore[[4]]
)

# Plot prediction model
mdlStudentOut$Att$predSlopesCore <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Att$lmeSlopesCore,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Att$predPltSlopesCore <-
    ggplot(data = mdlStudentOut$Att$predSlopesCore %>% filter(PID %in% studentOutPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesAttCore.png",
  mdlStudentOut$Att$predPltSlopesCore,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also find that the core need remains a strong predictor (even when letting the influence vary between participants).

TL;DR: The random slope adds significantly to the prediction model.

Need fulfillment and Interaction Quality

Based on the assertion that the relationship between core need fulfillment and outgroup attitudes is explained by a higher perceived interaction, the core need fulfillment should also significantly predict perceived interaction quality.

Unconditional model

Given that we now have the perceived interaction quality as our outcome variable of interest we again begin with an unconditional model (i.e., empty random intercept model), to see whether there is enough variance to explain within the participants. Similarly to before this is again done within the subsample of measurements during which an outgroup interaction was reported.

\[\begin{equation} \tag{24} InteractionQuality_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Qlt$lmerNull <-
  lme4::lmer(quality_overall ~ 1 + (1 | PID), 
             data = studentOutgroupInteraction) # use optim if it does not converge
mdlStudentOut$Qlt$lmeNull <-
  lme(
    quality_overall ~ 1,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.Qual.ML.r) #or with the lme function
summ(mdlStudentOut$Qlt$lmerNull, digits = 3, center = TRUE)
Observations 935
Dependent variable quality_overall
Type Mixed effects linear regression
AIC 8179.693
BIC 8194.215
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.144
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 78.849 1.017 77.536 100.065 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 7.506
Residual 18.318
Grouping Variables
Group # groups ICC
PID 108 0.144 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch~1 + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlStudentOut$Qlt$varNull <- 
  VarCorr(mdlStudentOut$Qlt$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlStudentOut$Qlt$tauNull <- 
  as.numeric(mdlStudentOut$Qlt$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlStudentOut$Qlt$sigmaNull <- 
  as.numeric(mdlStudentOut$Qlt$varNull[2])
# The ICC estimate (between/between+within) is:
mdlStudentOut$Qlt$IccNull <-
  (as.numeric(mdlStudentOut$Qlt$varNull[1]) / (as.numeric(mdlStudentOut$Qlt$varNull[1]) + as.numeric(mdlStudentOut$Qlt$varNull[2])))
mdlStudentOut$Qlt$IccPercNull <-
  ((as.numeric(mdlStudentOut$Qlt$varNull[1]) / (as.numeric(mdlStudentOut$Qlt$varNull[1]) + as.numeric(mdlStudentOut$Qlt$varNull[2])))) * 100

We again find a reasonable amount of variance within the participants.

Note to self: For completeness, we should probably check this against an empty model without random intercept.

random intercept with level one predictor

We again add the core interaction need fulfillment to the multilevel random intercept model.

\[\begin{equation} \tag{11} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Qlt$lmeInterceptCore <-
  lme(
    quality_overall ~ KeyNeedFullfillment,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlStudentOut$Qlt$lmerInterceptCore <- lmer(quality_overall ~ KeyNeedFullfillment + (1 | PID), 
       data = studentOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable quality_overall
Type Mixed effects linear regression
AIC 8081.658
BIC 8101.021
Pseudo-R² (fixed effects) 0.106
Pseudo-R² (total) 0.221
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 78.915 77.091 80.739 84.803 102.614 0.000
KeyNeedFullfillment 0.348 0.284 0.413 10.550 932.928 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 6.678
Residual 17.389
Grouping Variables
Group # groups ICC
PID 108 0.129 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(quality_overall ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Qlt$lmeNull, 
      mdlStudentOut$Qlt$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 19: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Qlt\)lmeNull 1 3 8180 8194 -4087
mdlStudentOut\(Qlt\)lmeInterceptCore 2 4 8082 8101 -4037 1 vs 2 100.035 < .001 
# Save variances
mdlStudentOut$Qlt$varInterceptCore <-
  lme4::VarCorr(mdlStudentOut$Qlt$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Qlt$varBtwInterceptCore <- 
  1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]) / as.numeric(mdlStudentOut$Qlt$varNull[1]))
mdlStudentOut$Qlt$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]) / as.numeric(mdlStudentOut$Qlt$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Qlt$varWithinInterceptCore <-
  1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]) / as.numeric(mdlStudentOut$Qlt$varNull[2]))
mdlStudentOut$Qlt$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]) / as.numeric(mdlStudentOut$Qlt$varNull[2]))) * 100

The predictor again adds a significant amount of explained variances beyond the empty model and looking at the slope coefficient, we find that the situational core need fulfillment relates significantly to perceived interaction quality. The core need explained substantial variance between participants (20.84%). The variance explained within participants is also medium (9.88%).

random slope

As before, we check whether further letting the effect of core need fulfillment vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{12} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudentOut$Qlt$lmeSlopesCore <-
  lme(
    quality_overall ~
      KeyNeedFullfillment,
    random = ~ 1 + KeyNeedFullfillment | PID,
    control = lmeControl(opt = "optim"),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Qlt$lmerSlopesCore <-
    lmer(
      quality_overall ~
        KeyNeedFullfillment +
        (1 + KeyNeedFullfillment | PID),
      data = studentOutgroupInteraction
    ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable quality_overall
Type Mixed effects linear regression
AIC 8047.231
BIC 8076.275
Pseudo-R² (fixed effects) 0.143
Pseudo-R² (total) 0.335
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 78.882 77.114 80.651 87.409 102.262 0.000
KeyNeedFullfillment 0.416 0.306 0.526 7.400 52.447 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 6.080
PID KeyNeedFullfillment 0.351
Residual 16.525
Grouping Variables
Group # groups ICC
PID 108 0.119 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(mdlStudentOut$Qlt$lmerSlopesCore,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Qlt$lmeInterceptCore, 
      mdlStudentOut$Qlt$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 20: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Qlt\)lmeInterceptCore 1 4 8082 8101 -4037
mdlStudentOut\(Qlt\)lmeSlopesCore 2 6 8047 8076 -4018 1 vs 2 38.42 < .001 
# Save variances
mdlStudentOut$Qlt$varSlopesCore <- 
  lme4::VarCorr(mdlStudentOut$Qlt$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Qlt$varBtwSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[1]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]))
mdlStudentOut$Qlt$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[1]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Qlt$varWithinSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[2]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]))
mdlStudentOut$Qlt$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[2]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]))) * 100

# Assumption Checks:
mdlStudentOut$Qlt$diagSlopesCore <-
  sjPlot::plot_model(mdlStudentOut$Qlt$lmerSlopesCore, type = "diag")
grid.arrange(
  mdlStudentOut$Qlt$diagSlopesCore[[1]],
  mdlStudentOut$Qlt$diagSlopesCore[[2]]$`PID`,
  mdlStudentOut$Qlt$diagSlopesCore[[3]],
  mdlStudentOut$Qlt$diagSlopesCore[[4]]
)

# Plot prediction model
mdlStudentOut$Qlt$predSlopesCore <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Qlt$lmeSlopesCore,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Qlt$predPltSlopesCore <-
    ggplot(data = mdlStudentOut$Qlt$predSlopesCore %>% filter(PID %in% studentOutPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesCore.png",
  mdlStudentOut$Qlt$predPltSlopesCore,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model.

Need fulfillment, Interaction Quality, and Attitudes

In our final main step, we will jointly consider the effect of core need fulfillment and perceived interaction quality on outgroup attitudes. We expect that if the relationship between core need fulfillment and outgroup attitudes is indeed explained by a higher perceived interaction quality, the interaction quality perception should assume the explained variance of the core contact need fulfillment.

random intercept with level one predictors

We thus add both the core need fulfillment and perceived interaction quality to a random intercept model of outgroup attitudes.

\[\begin{equation} \tag{13} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmeInterceptCoreQlt <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment + quality_overall,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    AttitudesDutch ~ KeyNeedFullfillment + quality_overall + (1 | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7179.308
BIC 7203.511
Pseudo-R² (fixed effects) 0.032
Pseudo-R² (total) 0.747
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.720 67.611 73.829 44.581 102.807 0.000
KeyNeedFullfillment 0.045 0.006 0.084 2.251 846.076 0.025
quality_overall 0.151 0.115 0.187 8.114 841.537 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.855
Residual 9.421
Grouping Variables
Group # groups ICC
PID 108 0.739 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(
  mdlStudentOut$Att$lmeNull, 
  mdlStudentOut$Att$lmeInterceptCoreQlt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 21: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeNull 1 3 7250 7265 -3622
mdlStudentOut\(Att\)lmeInterceptCoreQlt 2 5 7179 7203 -3585 1 vs 2 74.762 < .001 
# Save variances
mdlStudentOut$Att$varInterceptCoreQlt <-
  lme4::VarCorr(mdlStudentOut$Att$lmeInterceptCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwInterceptCoreQlt <- 
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))
mdlStudentOut$Att$varBtwPercInterceptCoreQlt <- 
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinInterceptCoreQlt <-
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))
mdlStudentOut$Att$varWithinPercInterceptCoreQlt <-
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))) * 100

Unsurprisingly, the model with the two predictors adds significantly beyond the empty unconditional model. However, more importantly, looking at the coefficients, we find that the effect of core need fulfillment indeed is indeed strongly reduced and the variance is explained by the perceived interaction quality. The variance explained in outgroup attitudes is of medium effect size (between: 2.32%, within: 9.42%).

random slope

We again check whether further letting the effects vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{14} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudentOut$Att$lmeSlopesCoreQlt <-
  lme(
    AttitudesDutch ~
      KeyNeedFullfillment + quality_overall,
    random = ~ 1 + KeyNeedFullfillment + quality_overall | PID,
    control = lmeControl(opt = "optim"),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Att$lmerSlopesCoreQlt <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment + quality_overall +
      (1 + KeyNeedFullfillment + quality_overall | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7137.672
BIC 7186.077
Pseudo-R² (fixed effects) 0.039
Pseudo-R² (total) 0.787
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.582 67.471 73.692 44.472 102.581 0.000
KeyNeedFullfillment 0.039 -0.004 0.082 1.789 24.796 0.086
quality_overall 0.173 0.114 0.233 5.723 50.585 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.844
PID KeyNeedFullfillment 0.055
PID quality_overall 0.192
Residual 8.733
Grouping Variables
Group # groups ICC
PID 108 0.767 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(
  mdlStudentOut$Att$lmeNull,
  mdlStudentOut$Att$lmeInterceptCoreQlt,
  mdlStudentOut$Att$lmeSlopesCoreQlt
) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 22: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeNull 1 3 7250 7265 -3622
mdlStudentOut\(Att\)lmeInterceptCoreQlt 2 5 7179 7203 -3585 1 vs 2 74.762 < .001
mdlStudentOut\(Att\)lmeSlopesCoreQlt 3 10 7138 7187 -3559 2 vs 3 51.111 < .001 
# Save variances
mdlStudentOut$Att$varSlopesCoreQlt <- 
  lme4::VarCorr(mdlStudentOut$Att$lmeSlopesCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwSlopesCoreQlt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]))
mdlStudentOut$Att$varBtwPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinSlopesCoreQlt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]))
mdlStudentOut$Att$varWithinPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]))) * 100

# Assumption Checks:
mdlStudentOut$Att$diagSlopesCoreQlt <- 
  sjPlot::plot_model(mdlStudentOut$Att$lmerSlopesCoreQlt, type = "diag")
grid.arrange(
  mdlStudentOut$Att$diagSlopesCoreQlt[[1]],
  mdlStudentOut$Att$diagSlopesCoreQlt[[2]]$`PID`,
  mdlStudentOut$Att$diagSlopesCoreQlt[[3]],
  mdlStudentOut$Att$diagSlopesCoreQlt[[4]]
)

# Plot prediction model
mdlStudentOut$Att$predSlopesCoreQlt <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Att$lmeSlopesCoreQlt,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Att$predPltSlopesCoreQlt <-
    ggplot(data = mdlStudentOut$Att$predSlopesCoreQlt %>% filter(PID %in% studentOutPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesAttCoreQlt.png",
  mdlStudentOut$Att$predPltSlopesCoreQlt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also find that the perceived interaction quality remains a strong predictor (even when letting the slopes vary between participants).

Check for robustness

To build further confidence in our results, we assess a few additional models that might offer alternative explanations of the effects we find.

Interaction Type

To make certain that the effect of core need fulfillment is specific to the interaction we compare the the effect to fulfillment of the situation core need when no outgroup interaction took place.

random intercept

Here we go back to the full dataset and add generalized situational core need fulfillment (either during an interaction or about the daytime in general) and whether an outgroup interaction happened as well as their interaction term to a random intercept model of outgroup attitudes.

\[\begin{equation} \tag{15} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudent$lmeInterceptAttCoreInt <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment * OutgroupInteraction,
    random =  ~ 1 | PID,
    data = studentInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlStudent$lmerInterceptAttCoreInt <- lmer(
    AttitudesDutch ~ KeyNeedFullfillment * OutgroupInteraction + (1 | PID),
    data = studentInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36668.990
BIC 36708.051
Pseudo-R² (fixed effects) 0.004
Pseudo-R² (total) 0.802
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.772 63.373 70.171 38.501 112.066 0.000
KeyNeedFullfillment 0.025 0.012 0.039 3.587 4859.664 0.000
OutgroupInteraction 2.638 1.915 3.360 7.154 4861.666 0.000
KeyNeedFullfillment:OutgroupInteraction 0.049 0.014 0.085 2.742 4853.890 0.006
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.349
Residual 9.152
Grouping Variables
Group # groups ICC
PID 113 0.801 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentInteractionType),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeAttNull, 
      mdlStudent$lmeInterceptAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 23: Student: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeAttNull 1 3 36741 36760 -18367
lmeInterceptAttCoreInt 2 6 36669 36708 -18328 1 vs 2 77.936 < .001 
# Save variances
mdlStudent$varInterceptAttCoreInt <- 
  lme4::VarCorr(mdlStudent$lmeInterceptAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwInterceptAttCoreInt <-
  1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[1]) / as.numeric(mdlStudent$varAttNull[1]))
mdlStudent$varBtwPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[1]) / as.numeric(mdlStudent$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinInterceptAttCoreInt <-
  1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[2]) / as.numeric(mdlStudent$varAttNull[2]))
mdlStudent$varWithinPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[2]) / as.numeric(mdlStudent$varAttNull[2]))) * 100

We find that the model explains significantly more variance than the empty null model. However, more interestingly, looking at the coefficients, we find that, as seen earlier, having an outgroup interaction has a strong effect on outgroup attitudes. Importantly, we find that there is a main effect of key need fulfillment but also a significant interaction effect of core need fulfillment and outgroup contact. This indicates that it is not simply key need fulfillment in general — but especially key need fulfillment during an outgroup contact that predicts more positive outgroup attitudes.

random slope

We again check whether further letting the effects vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{16} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudent$lmeSlopesAttCoreInt <- lme(
  AttitudesDutch ~
    KeyNeedFullfillment * OutgroupInteraction,
  random = ~ 1 + KeyNeedFullfillment + OutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = studentInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudent$lmerSlopesAttCoreInt <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment * OutgroupInteraction +
      (1 + KeyNeedFullfillment + OutgroupInteraction | PID),
    data = studentInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36473.506
BIC 36545.117
Pseudo-R² (fixed effects) 0.006
Pseudo-R² (total) 0.815
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.773 63.299 70.248 37.666 113.093 0.000
KeyNeedFullfillment 0.029 0.011 0.047 3.166 84.092 0.002
OutgroupInteraction 2.898 1.365 4.431 3.705 94.703 0.000
KeyNeedFullfillment:OutgroupInteraction 0.060 0.021 0.099 3.039 1342.328 0.002
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.752
PID KeyNeedFullfillment 0.055
PID OutgroupInteraction 6.932
Residual 8.794
Grouping Variables
Group # groups ICC
PID 113 0.820 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeAttNull, 
      mdlStudent$lmeInterceptAttCoreInt,
      mdlStudent$lmeSlopesAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  add_rownames(., var = "Description") %>%
  mutate(Description = gsub(".*\\$", "", Description)) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = TRUE,
    align = rep("l", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 24: Student: Model Comparison
Description Model df AIC BIC logLik Test L.Ratio p-value
lmeAttNull 1 3 36741 36760 -18367
lmeInterceptAttCoreInt 2 6 36669 36708 -18328 1 vs 2 77.936 < .001
lmeSlopesAttCoreInt 3 11 36474 36545 -18226 2 vs 3 205.459 < .001 
# Save variances
mdlStudent$varSlopesAttCoreInt <- 
  lme4::VarCorr(mdlStudent$lmeSlopesAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwSlopesAttCoreInt <- 
  1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[1]) / as.numeric(mdlStudent$varInterceptAttCoreInt[1]))
mdlStudent$varBtwPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[1]) / as.numeric(mdlStudent$varInterceptAttCoreInt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinSlopesAttCoreInt <- 
  1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[2]) / as.numeric(mdlStudent$varInterceptAttCoreInt[2]))
mdlStudent$varWithinPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[2]) / as.numeric(mdlStudent$varInterceptAttCoreInt[2]))) * 100

# Assumption Checks:
mdlStudent$diagSlopesAttCoreInt <-
  sjPlot::plot_model(mdlStudent$lmerSlopesAttCoreInt, type = "diag")
grid.arrange(
  mdlStudent$diagSlopesAttCoreInt[[1]],
  mdlStudent$diagSlopesAttCoreInt[[2]]$`PID`,
  mdlStudent$diagSlopesAttCoreInt[[3]],
  mdlStudent$diagSlopesAttCoreInt[[4]]
)

# Plot prediction model
mdlStudent$predSlopesAttCoreInt <- 
  studentInteractionType %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudent$lmeSlopesAttCoreInt,
                           studentInteractionType,
                           re.form = NA
                           )
         )

(
  mdlStudent$predPltSlopesAttCoreInt <-
    ggplot(data = mdlStudent$predSlopesAttCoreInt %>% filter(PID %in% studentPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Student_PredictionPlot_SlopesAttCoreInt.png",
  mdlStudent$predPltSlopesAttCoreInt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also see that when taking the possibility to varying slopes into account, the coefficient interpretations remains consistent (i.e., outgroup contact and its interaction with core need fulfillment remain important predictors of positive outgroup attitudes).

Plot Interaction

Before we move on, we shortly illustrate the interaction effect of how the effect of core need fulfillment differed by whether an outgroup contact took place or not. To this end we illustrate (1) the raw data points (without taking the nested nature into account), as well as a plot of the model predicted values and their prediction interval (taking the nested structure into account based; similar to an interaction plot).

# visualize interaction
## Without ML structure
ggplot(data = studentInteractionType,
       aes(x = KeyNeedFullfillment,
           y = AttitudesDutch,
           fill = OutgroupInteraction)) +
  #geom_point()+
  geom_smooth(method = 'lm',
              aes(linetype = OutgroupInteraction),
              color = "black") +
  #facet_wrap(~PID, ncol = 6)+
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  #scale_colour_manual(values=c("grey20", "black"), name="Intergroup Contact")+
  scale_y_continuous(
    limits = c(50, 100),
    breaks = seq(50, 100, by = 10),
    position = "left"
  ) +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(
    title = "Without ML stucture",
    x = "Fulfillment Core Need",
    y = "Outgroup Attitudes",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact"
  ) +
  theme_Publication() +
  theme(legend.position = "bottom",
        legend.key.size = unit(1, "cm"))

## With ML structure
# create parameters for prediction
datNew = data.frame(
  OutgroupInteraction = factor(rep(c("No", "Yes"), each = 21)),
  KeyNeedFullfillment = rep(seq(0, 100, 5), 2),
  PID = 0
)

# Predict values, clean up and calculate SE
PI <-
  merTools::predictInterval(
    merMod = mdlStudent$lmerSlopesAttCoreInt,
    newdata = datNew,
    level = 0.95,
    stat = "mean",
    type = "linear.prediction",
    include.resid.var = F,
    fix.intercept.variance = F
  )
mdlStudent$predInterceptAttCoreIntX <- 
  cbind(datNew, PI)
mdlStudent$predInterceptAttCoreIntX$se <-
  (mdlStudent$predInterceptAttCoreIntX$upr - mdlStudent$predInterceptAttCoreIntX$fit) / 1.96
rm(datNew, PI)

# Plot predicted values with SE
ggplot(
  mdlStudent$predInterceptAttCoreIntX,
  aes(x = KeyNeedFullfillment,
      y = fit,
      fill = OutgroupInteraction)
)+
  #geom_point() +
  geom_line(aes(linetype = as.factor(OutgroupInteraction)), size = 1) +
  #facet_wrap(~PID, ncol = 6)+
  geom_ribbon(data = mdlStudent$predInterceptAttCoreIntX,
              aes(ymin = fit - se, ymax = fit + se),
              alpha = 0.3) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  labs(
    x = "Fulfillment Core Need",
    y = "Outgroup Attitude (NL)",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact",
    title = "Based on Model Predictions"
  ) +
  theme_Publication()

# #### Bayesian estimation !! ONLY RUN ON FINAL RENDER !! Takes forever ####
# options(mc.cores = parallel::detectCores())  # Run many chains simultaneously
# brmfit <- brm(
#   AttitudesDutch ~ KeyNeedFullfillment * OutgroupInteraction + 
#     (1 + KeyNeedFullfillment + OutgroupInteraction | PID),
#   data = studentInteractionType,
#   family = gaussian,
#   iter = 1000,
#   chains = 4
# )
# 
# # create parameters for prediction:
# datNew = data.frame(
#   OutgroupInteraction = rep(c("No", "Yes"), each = 101),
#   KeyNeedFullfillment = rep(seq(0, 100, 1), 2)
# )
# # Save predicted values and adjust names and labels
# fitavg <-
#   cbind(datNew,
#         fitted(brmfit, newdata = datNew, re_formula = NA)[, -2])
# names(fitavg)[names(fitavg) == "Estimate"] = "pred"
# fitavg$se <- (fitavg$Q97.5 - fitavg$pred) / 1.96
# 
# # Plot Bayesian SE prediction interval
# ggplot(fitavg,
#        aes(x = KeyNeedFullfillment,
#            y = pred,
#            fill = OutgroupInteraction)) +
#   scale_x_continuous(breaks = seq(0, 100, 10)) +
#   scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
#   geom_line(aes(linetype = OutgroupInteraction), size = 1) +
#   geom_ribbon(aes(ymin = pred - se, ymax = pred + se), alpha = 0.2) +
#   scale_linetype_manual(values = c("dashed", "solid")) +
#   scale_fill_manual(values = c("darkgrey", "black")) +
#   labs(
#     x = "Fulfillment Core Need",
#     y = "Outgroup Attitude (NL)",
#     fill = "Intergroup Contact",
#     linetype = "Intergroup Contact",
#     title = "Based on Bayesian Prediction Interval"
#   ) +
#   theme_Publication()
# 
# # plot all overlayed posteriors:
# pst <- posterior_samples(brmfit, "b")
# ggplot(studentInteractionType,
#        aes(x = KeyNeedFullfillment, y = AttitudesDutch)) +
#   geom_point(shape = 4, alpha = .1) +
#   geom_tile() +
#   geom_abline(
#     data = pst,
#     aes(intercept = b_Intercept, slope = b_KeyNeedFullfillment),
#     alpha = .025,
#     size = .4
#   ) +
#   labs(title = "slope Posteriors",
#        x = "Fulfillment Core Need",
#        y = "Outgroup Attitudes (NL)") +
#   theme_Publication()
# rm(datNew, brmfit, fitavg, pst)

The plots indicate that especially once we take the nested data structure into account we can see a substantially stronger effect of core need fulfillment on outgroup attitudes during outgroup contacts than without outgroup contacts.

Other psychological needs

In a final step we check whether during the interaction the core situational need is a meaningful predictor even when taking other fundamental psychological needs into account. We focus on the three commonly considered self determination needs: competence, autonomy, and relatedness.

random intercept with level ome predictors

We add the core need fulfillment with the three self determination needs to a random intercept model of outgroup attitudes.

\[\begin{equation} \tag{17} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmeInterceptCoreSdt <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction,
    na.action = na.exclude
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlStudentOut$Att$lmerInterceptCoreSdt <- lmer(
    AttitudesDutch ~ KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction + (1 | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7212.060
BIC 7245.944
Pseudo-R² (fixed effects) 0.028
Pseudo-R² (total) 0.745
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.753 67.624 73.881 44.325 102.929 0.000
KeyNeedFullfillment 0.062 0.023 0.101 3.147 845.995 0.002
Competence 0.056 0.013 0.098 2.575 852.148 0.010
Autonomy 0.034 -0.015 0.083 1.344 870.215 0.179
RelatednessInteraction 0.058 0.033 0.084 4.425 848.412 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.943
Residual 9.519
Grouping Variables
Group # groups ICC
PID 108 0.737 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# To be compared against a model with only SDT needs
mdlStudentOut$Att$lmeInterceptSdt <-
  lme(
    AttitudesDutch ~ Competence + Autonomy + RelatednessInteraction,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction,
    na.action = na.exclude
  )

summ(
  mdlStudentOut$Att$lmerInterceptSdt <- lmer(
    AttitudesDutch ~ Competence + Autonomy + RelatednessInteraction + (1 | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7213.914
BIC 7242.958
Pseudo-R² (fixed effects) 0.025
Pseudo-R² (total) 0.744
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.702 67.558 73.847 44.072 102.829 0.000
Competence 0.067 0.025 0.109 3.122 853.478 0.002
Autonomy 0.034 -0.015 0.083 1.350 871.214 0.177
RelatednessInteraction 0.063 0.037 0.089 4.788 847.980 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 16.024
Residual 9.564
Grouping Variables
Group # groups ICC
PID 108 0.737 
# Compare new model to previous step
anova(
  mdlStudentOut$Att$lmeNull,
  mdlStudentOut$Att$lmeInterceptSdt, 
  mdlStudentOut$Att$lmeInterceptCoreSdt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 25: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeNull 1 3 7250 7265 -3622
mdlStudentOut\(Att\)lmeInterceptSdt 2 6 7214 7243 -3601 1 vs 2 42.156 < .001
mdlStudentOut\(Att\)lmeInterceptCoreSdt 3 7 7212 7246 -3599 2 vs 3 3.854 0.05 
# Save variances
mdlStudentOut$Att$varInterceptCoreSdt <-
  lme4::VarCorr(mdlStudentOut$Att$lmeInterceptCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwInterceptCoreSdt <- 
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))
mdlStudentOut$Att$varBtwPercInterceptCoreSdt <- 
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinInterceptCoreSdt <-
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))
mdlStudentOut$Att$varWithinPercInterceptCoreSdt <-
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))) * 100

We find that the the model with the added predictor indeed explains more variance in outgroup attitudes than the unconditional model and we find that adding the core need adds further explained variance — beyond the self determination needs. Looking at the individual coefficients, we find that the situational core need relates significantly to outgroup attitudes, that it is a stronger predictor than any of the self determination theory needs and that it assumes some of the variance explained by the self determination theory needs (when compared to a model without the core need).

random slope

We again check whether further letting the effects vary between participants explains additional variance in outgroup attitudes (i.e., random slope).

\[\begin{equation} \tag{18} Attitude_{ti} = \\gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Check Correlations 
studentOutgroupInteraction %>% 
  select(AttitudesDutch, KeyNeedFullfillment, Competence, Autonomy, RelatednessInteraction) %>%
  pairs.panels.new

studentRed <- 
  studentOutgroupInteraction %>% 
  # select(PID, AttitudesDutch, KeyNeedFullfillment, Competence, Autonomy, RelatednessInteraction) %>%
  group_by(PID) %>%
  filter(
    sd(AttitudesDutch) != 0,
    sd(KeyNeedFullfillment) != 0,
    sd(Competence) != 0,
    sd(Autonomy) != 0,
    sd(RelatednessInteraction) != 0
  ) %>%
  ungroup

studentRedCor <- 
  studentOutgroupInteraction %>%
  group_by(PID) %>%
  summarise(
    rAttCore = cor(AttitudesDutch, KeyNeedFullfillment),
    rAttComp = cor(AttitudesDutch, Competence),
    rAttAut = cor(AttitudesDutch, Autonomy),
    rAttRel = cor(AttitudesDutch, RelatednessInteraction),
    rCoreComp = cor(KeyNeedFullfillment, Competence),
    rCoreAut = cor(KeyNeedFullfillment, Autonomy),
    rCoreRel = cor(KeyNeedFullfillment, RelatednessInteraction),
    rCompAut = cor(Competence, Autonomy),
    rCompRel = cor(Competence, RelatednessInteraction),
    rAutRel = cor(Autonomy, RelatednessInteraction)
  ) %>%
  filter_at(vars(-PID), all_vars(abs(.) != "1"))
  # mutate(corMean = rowMeans(abs(.[2:ncol(.)]))) %>%
  # filter(corMean != "1")

studentRed2 <- 
  studentOutgroupInteraction %>%
  filter(PID %in% studentRedCor$PID)

# Create and save Model (optimizer needed to reach convergence) 
# Problem because some PPTs have 100 on all measures (SD = 0) AND/OR
# For some all cor = 1 or -1
# Removing these PPTs or their measurement beeps doesn't help
# However, removing eithe the Core Need, Autonomy, or Relatedness from the random slopes lets the model converge
# --> looks like there is an overfitting issue? But why does it work in lmer?
# FOR NOW: Autonomy removed from random slopes because weakest predictor
mdlStudentOut$Att$lmeSlopesCoreSdt <-
  nlme::lme(
    AttitudesDutch ~
      KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction,
    random = ~ 1 + KeyNeedFullfillment + Competence + RelatednessInteraction | PID, # Autonomy + 
    #method = "ML",
    #control = lmeControl(opt = "optim", optimMethod = "L-BFGS-B"),
    control = lmeControl(opt = "optim", maxIter = 100, msMaxIter = 100),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Att$lmerSlopesCoreSdt <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction +
      (1 + KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7214.338
BIC 7315.990
Pseudo-R² (fixed effects) 0.034
Pseudo-R² (total) 0.765
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.672 67.578 73.766 44.775 103.278 0.000
KeyNeedFullfillment 0.085 0.031 0.139 3.072 40.508 0.004
Competence 0.052 0.008 0.096 2.303 92.107 0.024
Autonomy 0.027 -0.022 0.076 1.083 254.850 0.280
RelatednessInteraction 0.066 0.035 0.097 4.180 40.814 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.703
PID KeyNeedFullfillment 0.124
PID Competence 0.047
PID Autonomy 0.022
PID RelatednessInteraction 0.069
Residual 9.089
Grouping Variables
Group # groups ICC
PID 108 0.749 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Att$lmeInterceptCoreSdt, 
      mdlStudentOut$Att$lmeSlopesCoreSdt) %>%
  as.data.frame() %>%
  select(-call) %>%
  mutate(
    L.Ratio = round(L.Ratio, 3),
    `p-value` = ifelse(`p-value`>=.001, round(`p-value`, 3), "< .001")
  ) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 26: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeInterceptCoreSdt 1 7 7212 7246 -3599
mdlStudentOut\(Att\)lmeSlopesCoreSdt 2 16 7207 7284 -3587 1 vs 2 23.1 0.006 
# Save variances
mdlStudentOut$Att$varSlopesCoreSdt <- 
  lme4::VarCorr(mdlStudentOut$Att$lmeSlopesCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwSlopesCoreSdt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]))
mdlStudentOut$Att$varBtwPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinSlopesCoreSdt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]))
mdlStudentOut$Att$varWithinPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]))) * 100

# Assumption Checks:
mdlStudentOut$Att$diagSlopesCoreSdt <- 
  sjPlot::plot_model(mdlStudentOut$Att$lmerSlopesCoreSdt, type = "diag")
grid.arrange(
  mdlStudentOut$Att$diagSlopesCoreSdt[[1]],
  mdlStudentOut$Att$diagSlopesCoreSdt[[2]]$`PID`,
  mdlStudentOut$Att$diagSlopesCoreSdt[[3]],
  mdlStudentOut$Att$diagSlopesCoreSdt[[4]]
)

# Plot prediction model
mdlStudentOut$Att$predSlopesCoreSdt <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Att$lmeSlopesCoreSdt,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Att$predPltSlopesCoreSdt <-
    ggplot(data = mdlStudentOut$Att$predSlopesCoreSdt %>% filter(PID %in% studentOutPltIDs),
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesAttCoreStd.png",
  mdlStudentOut$Att$predPltSlopesCoreSdt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

We find that adding the random slopes does add significantly beyond the random intercept model. We also see that when taking the possibility to varying slopes into account, the coefficient interpretations remains consistent (i.e., core need and relatedness remain the strongest and only significant predictors).

Young Medical Professional Sample

Contact Hypothesis

Allport’s Conditions

Need Fulfillment

Software Information

The full session information with all relevant system information and all loaded and installed packages is available in the collapsible section below.

System Info
Table 27: R environment session info for reproducibility of results
Setting Value
version R version 4.1.1 (2021-08-10)
os macOS Big Sur 10.16
system x86_64, darwin17.0
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/Amsterdam
date 2021-10-20

Package Info
Table 28: Package info for reproducibility of results
Package Loaded version Date Source
bookdown 0.24 2021-09-02 CRAN (R 4.1.0)
brms 2.16.1 2021-08-23 CRAN (R 4.1.0)
data.table 1.14.0 2021-02-21 CRAN (R 4.1.0)
devtools 2.4.2 2021-06-07 CRAN (R 4.1.0)
dplyr 1.0.7 2021-06-18 CRAN (R 4.1.0)
ellipse 0.4.2 2020-05-27 CRAN (R 4.1.0)
Formula 1.2-4 2020-10-16 CRAN (R 4.1.0)
ggpattern 0.2.0 2021-10-11 Github ()
ggplot2 3.3.5 2021-06-25 CRAN (R 4.1.0)
ggthemes 4.2.4 2021-01-20 CRAN (R 4.1.0)
gridExtra 2.3 2017-09-09 CRAN (R 4.1.0)
gtsummary 1.4.2 2021-07-13 CRAN (R 4.1.0)
haven 2.4.3 2021-08-04 CRAN (R 4.1.0)
Hmisc 4.5-0 2021-02-28 CRAN (R 4.1.0)
jtools 2.1.4 2021-09-03 CRAN (R 4.1.0)
kableExtra 1.3.4 2021-02-20 CRAN (R 4.1.0)
knitr 1.36 2021-09-29 CRAN (R 4.1.0)
lattice 0.20-44 2021-05-02 CRAN (R 4.1.1)
lme4 1.1-27.1 2021-06-22 CRAN (R 4.1.0)
lubridate 1.7.10 2021-02-26 CRAN (R 4.1.0)
mada 0.5.10 2020-05-25 CRAN (R 4.1.0)
Matrix 1.3-4 2021-06-01 CRAN (R 4.1.1)
mvmeta 1.0.3 2019-12-10 CRAN (R 4.1.0)
mvtnorm 1.1-2 2021-06-07 CRAN (R 4.1.0)
nlme 3.1-152 2021-02-04 CRAN (R 4.1.1)
pander 0.6.4 2021-06-13 CRAN (R 4.1.0)
papaja 0.1.0.9997 2021-10-11 Github ()
plotly 4.9.4.9000 2021-08-28 Github ()
plyr 1.8.6 2020-03-03 CRAN (R 4.1.0)
psych 2.1.6 2021-06-18 CRAN (R 4.1.0)
RColorBrewer 1.1-2 2014-12-07 CRAN (R 4.1.0)
Rcpp 1.0.7 2021-07-07 CRAN (R 4.1.0)
remedy 0.1.0 2018-12-03 CRAN (R 4.1.0)
reshape2 1.4.4 2020-04-09 CRAN (R 4.1.0)
rmarkdown 2.11 2021-09-14 CRAN (R 4.1.1)
sessioninfo 1.1.1 2018-11-05 CRAN (R 4.1.0)
stringi 1.7.5 2021-10-04 CRAN (R 4.1.1)
stringr 1.4.0 2019-02-10 CRAN (R 4.1.0)
survival 3.2-12 2021-08-13 CRAN (R 4.1.1)
tibble 3.1.5 2021-09-30 CRAN (R 4.1.0)
tidyr 1.1.4 2021-09-27 CRAN (R 4.1.0)
usethis 2.0.1 2021-02-10 CRAN (R 4.1.0)

Full Session Info (including loaded but unattached packages — for troubleshooting only)

R version 4.1.1 (2021-08-10)

Platform: x86_64-apple-darwin17.0 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages:

  • grid
  • stats
  • graphics
  • grDevices
  • datasets
  • utils
  • methods
  • base

other attached packages:

  • lubridate(v.1.7.10)
  • reshape2(v.1.4.4)
  • stringi(v.1.7.5)
  • stringr(v.1.4.0)
  • papaja(v.0.1.0.9997)
  • kableExtra(v.1.3.4)
  • Hmisc(v.4.5-0)
  • Formula(v.1.2-4)
  • survival(v.3.2-12)
  • lattice(v.0.20-44)
  • tidyr(v.1.1.4)
  • dplyr(v.1.0.7)
  • plyr(v.1.8.6)
  • data.table(v.1.14.0)
  • mada(v.0.5.10)
  • mvmeta(v.1.0.3)
  • ellipse(v.0.4.2)
  • mvtnorm(v.1.1-2)
  • devtools(v.2.4.2)
  • usethis(v.2.0.1)
  • pander(v.0.6.4)
  • tibble(v.3.1.5)
  • sessioninfo(v.1.1.1)
  • gtsummary(v.1.4.2)
  • jtools(v.2.1.4)
  • nlme(v.3.1-152)
  • lme4(v.1.1-27.1)
  • Matrix(v.1.3-4)
  • ggpattern(v.0.2.0)
  • gridExtra(v.2.3)
  • plotly(v.4.9.4.9000)
  • RColorBrewer(v.1.1-2)
  • haven(v.2.4.3)
  • ggthemes(v.4.2.4)
  • ggplot2(v.3.3.5)
  • psych(v.2.1.6)
  • brms(v.2.16.1)
  • Rcpp(v.1.0.7)
  • bookdown(v.0.24)
  • remedy(v.0.1.0)
  • knitr(v.1.36)
  • rmarkdown(v.2.11)

loaded via a namespace (and not attached):

  • estimability(v.1.3)
  • coda(v.0.19-4)
  • dygraphs(v.1.1.1.6)
  • rpart(v.4.1-15)
  • inline(v.0.3.19)
  • generics(v.0.1.0)
  • callr(v.3.7.0)
  • commonmark(v.1.7)
  • proxy(v.0.4-26)
  • chron(v.2.3-56)
  • tzdb(v.0.1.2)
  • webshot(v.0.5.2)
  • xml2(v.1.3.2)
  • httpuv(v.1.6.3)
  • StanHeaders(v.2.21.0-7)
  • assertthat(v.0.2.1)
  • xfun(v.0.26)
  • hms(v.1.1.1)
  • jquerylib(v.0.1.4)
  • bayesplot(v.1.8.1)
  • evaluate(v.0.14)
  • promises(v.1.2.0.1)
  • fansi(v.0.5.0)
  • igraph(v.1.2.6)
  • DBI(v.1.1.1)
  • tmvnsim(v.1.0-2)
  • htmlwidgets(v.1.5.4)
  • tensorA(v.0.36.2)
  • stats4(v.4.1.1)
  • purrr(v.0.3.4)
  • ellipsis(v.0.3.2)
  • crosstalk(v.1.1.1)
  • backports(v.1.2.1)
  • V8(v.3.4.2)
  • insight(v.0.14.3)
  • markdown(v.1.1)
  • RcppParallel(v.5.1.4)
  • vctrs(v.0.3.8)
  • remotes(v.2.4.0)
  • sjlabelled(v.1.1.8)
  • abind(v.1.4-5)
  • cachem(v.1.0.6)
  • withr(v.2.4.2)
  • checkmate(v.2.0.0)
  • emmeans(v.1.6.3)
  • xts(v.0.12.1)
  • prettyunits(v.1.1.1)
  • mnormt(v.2.0.2)
  • svglite(v.2.0.0)
  • cluster(v.2.1.2)
  • lazyeval(v.0.2.2)
  • crayon(v.1.4.1)
  • labeling(v.0.4.2)
  • pkgconfig(v.2.0.3)
  • pkgload(v.1.2.1)
  • blme(v.1.0-5)
  • nnet(v.7.3-16)
  • rlang(v.0.4.11)
  • lifecycle(v.1.0.1)
  • miniUI(v.0.1.1.1)
  • colourpicker(v.1.1.0)
  • modelr(v.0.1.8)
  • distributional(v.0.2.2)
  • rprojroot(v.2.0.2)
  • matrixStats(v.0.60.1)
  • datawizard(v.0.2.0)
  • loo(v.2.4.1)
  • boot(v.1.3-28)
  • zoo(v.1.8-9)
  • base64enc(v.0.1-3)
  • gamm4(v.0.2-6)
  • ggridges(v.0.5.3)
  • processx(v.3.5.2)
  • png(v.0.1-7)
  • viridisLite(v.0.4.0)
  • parameters(v.0.14.0)
  • rootSolve(v.1.8.2.2)
  • arm(v.1.11-2)
  • readr(v.2.0.2)
  • jpeg(v.0.1-9)
  • shinystan(v.2.5.0)
  • ggeffects(v.1.1.1)
  • scales(v.1.1.1)
  • memoise(v.2.0.0)
  • magrittr(v.2.0.1)
  • threejs(v.0.3.3)
  • compiler(v.4.1.1)
  • rstantools(v.2.1.1)
  • cli(v.3.0.1)
  • lmerTest(v.3.1-3)
  • interactions(v.1.1.5)
  • TMB(v.1.7.21)
  • ps(v.1.6.0)
  • Brobdingnag(v.1.2-6)
  • htmlTable(v.2.2.1)
  • MASS(v.7.3-54)
  • mgcv(v.1.8-36)
  • tidyselect(v.1.1.1)
  • forcats(v.0.5.1)
  • mixmeta(v.1.1.3)
  • projpred(v.2.0.2)
  • highr(v.0.9)
  • yaml(v.2.2.1)
  • latticeExtra(v.0.6-29)
  • bridgesampling(v.1.1-2)
  • sass(v.0.4.0)
  • tools(v.4.1.1)
  • lmom(v.2.8)
  • parallel(v.4.1.1)
  • rstudioapi(v.0.13)
  • foreach(v.1.5.1)
  • foreign(v.0.8-81)
  • gld(v.2.6.2)
  • posterior(v.1.1.0)
  • farver(v.2.1.0)
  • sjPlot(v.2.8.9)
  • digest(v.0.6.28)
  • shiny(v.1.6.0)
  • broom(v.0.7.9.9001)
  • performance(v.0.7.3)
  • later(v.1.3.0)
  • httr(v.1.4.2)
  • rsconnect(v.0.8.24)
  • effectsize(v.0.4.5)
  • sjstats(v.0.18.1)
  • colorspace(v.2.0-2)
  • rvest(v.1.0.1)
  • fs(v.1.5.0)
  • splines(v.4.1.1)
  • rematch2(v.2.1.2)
  • expm(v.0.999-6)
  • Exact(v.3.0)
  • renv(v.0.14.0)
  • shinythemes(v.1.2.0)
  • systemfonts(v.1.0.2)
  • xtable(v.1.8-4)
  • jsonlite(v.1.7.2)
  • nloptr(v.1.2.2.2)
  • rstan(v.2.21.2)
  • testthat(v.3.0.4)
  • gt(v.0.3.1)
  • R6(v.2.5.1)
  • broom.mixed(v.0.2.7)
  • pillar(v.1.6.3)
  • htmltools(v.0.5.2)
  • mime(v.0.12)
  • glue(v.1.4.2)
  • fastmap(v.1.1.0)
  • minqa(v.1.2.4)
  • DT(v.0.19)
  • class(v.7.3-19)
  • codetools(v.0.2-18)
  • pkgbuild(v.1.2.0)
  • utf8(v.1.2.2)
  • bslib(v.0.3.0)
  • numDeriv(v.2016.8-1.1)
  • pbkrtest(v.0.5.1)
  • curl(v.4.3.2)
  • DescTools(v.0.99.43)
  • gtools(v.3.9.2)
  • shinyjs(v.2.0.0)
  • glmmTMB(v.1.1.2.3)
  • merTools(v.0.5.2)
  • desc(v.1.3.0)
  • munsell(v.0.5.0)
  • e1071(v.1.7-9)
  • iterators(v.1.0.13)
  • broom.helpers(v.1.4.0)
  • sjmisc(v.2.8.7)
  • gtable(v.0.3.0)
  • bayestestR(v.0.10.5)


References